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ABSTRACT 


The main goal of this research is to model a flexible missile with structural] 
flexibility utilizing the Equivalent Rigid Link System (ERLS) with an enhanced 
natural mode discretization. Dynamic analysis of the flexible missile in 
2-Dimension motion is presented. 

Computer simulation is performed where the pitch angle of the missile is 
controlled with a rigid-body controller. The effects of increasing payloads and 


speed to the system performance are discussed. 


THESIS DISCLAIMER 


The reader is cautioned that computer programs developed in this research 
may not have been exercised for all cases of interest. While every effort has been 
made, within the time available, to ensure that the programs are free of compu- 
tational and logic errors, they cannot be considered validated. Any application 


of these programs without additional verification is at the risk of the user. 
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I. INTRODUCTION 


A. BACKGROUND 

With the introduction of long slender missiles such as the Vanguard, the 
Redstone, and various ballistic missiles, the problem of the structural flexibility 
became severe [Ref. 1]. Due to the limited thrust available from rocket engines, 
these missiles had to be as light as possible. This meant a sacrifice in structural 
rigidity. 

Missile flexure causes additional aerodynamic loads which in turn cause ad- 
ditional flexure. Coupling occurs between the elastic modes and the control sys- 
tem as the control system gyros sense the flexure motion and the rigid body 
motion. It has become necessary to actively control the flexible structure and 
thereby reduce the structural loads and improve the vehicle response such as po- 
sition, velocity and acceleration. Reduced structural loads will also offer potential 
for reduced bending stress and fatigue problems. 

The objective of this thesis is to develop a dynamic model for a flexible missile 
and study the dynamic behavior of the flexible missile. Simulation is a valuable 
tool in the design of new missile systems and in the modification or evaluation 
of existing systems. A missile simulation allows the engineer to evaluate his design 
without the expense of actually building and flying the actual missile. System 
dynamics can be investigated through simulation with a substantial savings in 


time and expense [Ref. 2]. 


B. LITERATURE REVIEW 

Flexible missile modeling centers on the relationship between the large, rigid- 
body motion and the small motions due to structural flexibility. 

Jenkins [Ref. 2] expresses some techniques used in deriving the equations of 
motion of a rigid missile for a six degree-of-freedom (6-DOF) simulation. The 
rigid missiles are characterized by their larger size and low ratio of payload to 
total weight. Rigid missile dynamic equations were developed using the 


Newton-Euler Method. The moments and forces along with the mass and mo- 


ments of inertia are assumed to be known in the body coordinate frame. The 
transformation between global’ and local coordinate frame is achieved with a 
non-homogenous coordinate transformation matrix [Ref. 3: pp. 342]. The result- 
ing rigid-body equations of motion produced the understanding for the derivation 
of the dynamic equations of the flexible missile. 

The Equivalent Rigid Link System [Ref. 4] describes the large-motion 
kinematics of the system by the ERLS motion and the small motion relative to 
the ERLS. The application of the finite element techniques and Lagrangian dy- 
namics produces two sets of coupled, nonlinear, ordinary differential equations 
of motion, of which one set is for the large motions and the other set for the small 
motions. The small motion is described by the superposition of vibration modes. 
The modes of the vibration of the flexible bar was derived with the simple-beam 
theory [Ref. 5: pp.221]. In simple beam theory, it is assumed that the rotation of 
the element is insignificant compared to the vertical translation and the shear 
deformation is small relative to the bending deformation. This assumption 1s valid 
if the ratio between the length of the bar and its height is relatively large (more 
than 10). The set of large motion equations are non-linear in both the large and 
Small motion variables while the set of small motion equations are linear in the 
small motion variable and non-linear in the large motion variables. A solution 
technique called the Sequential Integration Method [Ref. 6] was developed which 
allows efficient simulations of systems with inertia coupled motions having non- 
linear slow motion (large motion) with linear fast motion (small motion). The 
ERLS model presents a complete, efficient dynamic model able to describe large 
motion, smal! motion and their coupling. 

An ERLS model of a flexible spacecraft boom was developed and a computer 
Simulation was performed [Ref. 7]. The equations of motion were solved using 
the Sequential Integration Method. A spatial finite discretization of the boom 
Structure and the application of an assumed polynomial modal response were 
utilized in the approximate solution to the equations of motion. This work was 


the basis work for the modeling of the flexible missiles. 


Ganon [Ref. 8] performed an experimental validation of the ERLS dynamic 
model in a vertical plane motion. The small motion was modeled using a shape 
matrix derived from superposition of natural modes. The agreement between the 
simulation results and the experimental results justify the application of the 
ERLS using a natural-mode shape matrix to model the dynamic response of 


flexible missile. 


C. OUTLINE 

In this study, the ERLS dynamic model is used to derive the system equations 
of motion for a flexible missile in 2-D motion. Dynamic response is predicted by 
solving the equations of motion using the Sequential Integration Method. The 
application of an assumed natural mode shape and the spatial finite element 
discretization of the missile provide an approximate solution. Computer simu- 
lation for the flexible missile is performed using MATLAB on the MACINTOSH 
computer. A rigid body controller is included in the simulation to control the 
pitch angle of the flexible missile. 

The ERLS dynamic model of the flexible missile is presented in Chapter Two. 
A rigid body controller for the flexible missile is described in Chapter Three. The 
computer simulation methodology and simulation results are presented in Chap- 
ter Four. The conclusions and recommendations for future work are presented 


in Chapter Five. 


Il. MODEL FORMULATION FOR A FLEXIBLE MISSILE 


A. KINEMATICS 
In our model, 2-D inertial reference frame, i.e., X and Y, is used for the global 
motion as shown in Fig.l. The body-fixed coordinate frame, i.e., X and y, is se- 


lected to describe the missile local motion. 


x=Local xX axis 
y=Local y axis 
X=Global X axis 
Y=Global Y axis 
L=Missile Length 
D=Missile Diameter 





Figure 1. The general configuration of the missile 


The following assumptions were made for the flexible missile: 

(1) Material density is constant throughout the body and the steel was chosen 
for modeling. | 

(2) A uniform circular cross section is assumed. 

The geometric configuration parameters and material properties of the flexi- 


ble missile are listed as 


Diameter =0.12 meter 

Length =4.0 meter 

Material Density = 7861.05 kg/m? 

Young’s Modulus =2.0 x 10!' pascal 

The concept of the ERLS is applied to model the kinematics of flexible mis- 
sile. The main idea of the Equivalent Rigid Link System (Fig.2) is to separate the 


kinematics of the flexible body into large and small motions. 








| v(L 
| <a) 
v(0) = Missile Base Deflection 
(0) = Missile Base Slope 
v(L) = Missile Tip Deflection 
@(L) = Missile Tip Slope 
R = Global Position Vector 
. of Missile Base 
R, = Absolute Position Vector 
of the base of ERLS 
T = Local Position Vector 
d = Deformation Vector 
= Equivalent Rigid Link 
mz 3 = Theoretical Missile 





Position 
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' Figure 2. Equivalent Rigid Link System (ERLS) 


The large motion of the missile will be represented by missile base position in 
X direction Xj, missile base position in Y direction Y,, and missile pitch angle @. 
The small motions resulting from flexible motion are measured relative to the lo- 
cal coordinate frame x, y. v(O) and (0) are the nodal displacement and slope 
of the missile base respectively. The absolute (global) position of a point on the 
flexible missile is obtained from coordinate transformations. The global position 
(R) of any point position can be defined using ERLS in terms of a local position | 
vector (r), a deformation vector (d), and a coordinate transformation matrix (W), 


ree 
R=W +a) (AT 


The transformation matrix (W) between the large motion and small motion co- 


ordinate 1S 


1 0 0 
WW =| X, cos(@) sin(@) (2 =) 
Yo sin(@) cos(@) 


and the local rigid body position vector 1s 


I 


r= 


=| x (2 4 3) 
0 
The deformation vector (d) is expressed in terms of a nodal displacement 


vector U and a shape function @ as 


d=0U (2 — 4) 
where 
U=[v0) o(0)]" (2 — 5) 


The shape function @ was derived utilizing a natural-mode superposition and 
a finite element concept. The Finite Element Method (FEM) was utilized to 
discretize the flexible body displacements and assigning the nodes. Displacements 
are for each point along the missile, a function of location and time, and it is 
necessary to discretize the deformations to obtain an ordinary differential 
equation. The natural mode shape function of a beam is used to represent the 
flexural motion of the flexible missile. Only the first two mode shapes are used. 
The flexible missile is modeled as a continuous Euler-Bernoulli free-free beam, 
neglecting shear deformation and rotary inertia effects. The nodal points are the 
base points of the flexible missile. Appendix A shows the mode shape function 


derivation. In this case, ¢ is found as, 


0 0 

o—=!0 0 (2 — 6) 
a b 

Where a and b are defined as following, 

a= F,(C,(cos B)x + cosh £)x) + (sin Byx + sinh £,x)) 

+ F3(C,( cos B,x + cosh B,x) + (sin Bx + sinh f,x)) (2 —7) 

b = F,(C,( cos B,x + cosh £B,x) + (sin B)x + sinh £,-x)) 

+ F4(C,( cos B,x + cosh B,x) + (sin B,x + sinh £,x)) (2 — 8) 


Lagrangian dynamics is used to acquire the equations of motion and the 
kinetic energy of the system will be needed in the development of the Lagrangian 


expression. The absolute velocity is listed as follows, 


R=W(r +d)+ Wd | (2 — 9) 


B. KINETICS 
Lagrangian Dynamics is a systematic way to derive equations of motion for 
complex systems like flexible missiles. Lagrange’s equation is written as, 


d(0KE)_ @KE , aPE 


dt ( aq; = aq; éq; aa QO; (i i ee (2 oa 10) 


where 
K.E. = Kinetic Energy 
PE. = Potentiaeenersy, 
g,; = Generalized coordinates 


O, = Generalized forces 


a 


n = Number of degrees of freedom 


The generalized coordinates are chosen to be, 
q =[Xy Yo @ v(0) &(0)]" (2-11) 


The kinetic energy of the system is defined as follows, 


Ki s |RTRam QO. 


By substituting Eq.(2-1) into Eq. (2-12), the Kinetic energy is written as 


KE = > [WT Wr + 2r WT WOU + 2r'W WOU 


+ U6 WwW WoU +2U' 6 WwW WoU + U6 WT WoU dm (Que 


The potential energy of the system includes the strain energy of the flexible 


missile and the gravitational potential energy. 


P.E. = > [UT ’CrUdx — [R zd (2 — 14) 


Scieke 
[ 
C = Rigidity Matrix 
g = Gravitational acceleration vector 


Spatial derivative of the shape function 


Appendix B includes the development of the equations of motions using 
Lagrange equations. 

Generalized forces will be found by virtual work principle. It is assumed that 
_ the only force is the thrust force which is applied to the base of the flexible missile 
as shown in Fig.3. Other forces like aerodynamic and damping forces are neg- 
Iected. 


Equivalent Rigid Link. 
system 

Theoretical Missile 
Position 


Applied Force 





Figure 3. Applied Forces on Flexible Missile 


Ess force is composed of two components i.e., 


T,. = Tcos(6 + &(0)) (2 — 15) 


T, = Tsin(5 + (0) (2 — 16) 
and a moment applied on the missile 1s, 
Mp = —T cos(o + (0))v(0) (2 — 17) 


By applying virtual work principle, the generalized forces of flexible missile is 


found below. 


T cos(@) — Td sin(@) — T@(0) sin(@) 
F, =| Tsin(6) + Té cos(0) +T®(0) cos(6) (2 — 18) 
—Tv(0) 


= Td + TPO) 
F,= (2 — 19) 
—7Tv(0) 
where 
T = Force 


6 = Control Angle being assumed small 


Large motion generalized force vector 


ah 


Small motion generalized force vector 

Appendix C shows the derivation of the generalized forces using the Virtual 
Work Principle. 

The derivation of the equations of motion from the Lagrangian formulation 
yields two sets of coupled equations. One set describes the large motions and the 
another set describes the small motions. These two sets of equations are non- 


linear, coupled, second-order, ordinary differential equations of the form, 


My & + My, U = F, (2 — 20) 


—* 


Mrg & + Mn U+ G,U+K,U = F, O21) 


where: 

Mj, = 3x3 effective inertia matrix for large motions 

M,, = 3x2 coupled inertia matrix of the small motion effect on large 
motions 

F, = 3x1 load vector for the large motions 

M,qg = 2X3 coupled inertia matrix of the large motion effect on small 
motions 

Min = 2X2 effective inertia matrix for small motions 

G, = 2x2 > gyroscopic matrix 

K, = 2x2 stiffness matrix 

F- = 2x1 load vector for the small motions 

é = 3x1 vector, generalized coordinates of the large motions 

U = 2x1 vector, generalized coordinates of the small motions 


The detailed development of the equations of motion and definitions of the terms 


are described in Appendix B. 
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IW. THE DEVELOPMENT OF A RIGID-BODY CONTROLLER 


In our research, the controller is developed based on the rigid-body assump- 
tion. The purpose of including the rigid-body controller is to perform computer 
simulation and study the dynamic behavior of the flexible missile. The pitch an- 
gle of the missile is controlled by the rigid-body controller which is a single input 
Single output (SISO) controller. A general configuration of the flexible missile 


and rigid-body controller are shown in Fig.4. 


Control 
Rigid-Body} Angle 
Controller Missile 





Figure 4. The flexible missile with rigid-body controller 


The pitch angle is the control variable and the desired states are desired pitch 
angle @,, desired angular velocity @,, and desired angular acceleration @,. The 
control angle (6) is assumed small. 


The desired error function is defined as 
E+ K+ K,e=0 (3-1) 


where, 


K, = The position feedback gain. 
K, = The velocity feedback gain. 


€=0-—0, = The error in the position. 
é= 6 — 6, = The error in velocity. 
é=6-6 


@, = The error in acceleration. 
The control design begins with the equation of rigid-body (large motion) that 
is obtained by modifying and expanding Eq.(2-20) as 


M,, OM i h b 
jy al Lal Lalla 3 
My, Mp2} M2 ty hy b, 


Where, 
Mem) M12 
My) _— aa aa (3 _— 3) 
eo) Mq(2,2) 
Viale) 
Mpn=|.™ (3 — 4) 
Mo(2,3) 
M>) = ees!) Moa\2;2) 1] eC — 5) 
voy = ze) Al (3 — 6) 
ex 
a | | (3 — 7) 
Yo 
n2 = [6] (3 — 8) 


fi =Gravity, centrifugal and coriolis forces. 


fj; =Gravity, centrifugal and coriolis forces. 


h _ T cos(@) 
7 T sin(@) 


hy =(0} 
- ae 
b, == 

T cos(@) 
by = [0] 


Eq.(3-2) can be rewritten in a tensor form as 
My) + Min =f, + hy + 6,6 


Mn) + Mon, =f, 
Prom: Ea.(3-1 1), 
1 = Mi fi + hy + 5,0 = Miia| 


Substituting Eq.(3-13) into Eq.(3-12), we find, 


1 
A= M2, — MyM\, My 

Dee ee ae 
fp ean lf 1 hy | 


<b = pe 
B — M>,M), b, 


(3 ee 
(3 — 10) 
(3-11) 
(3 — 12) 
(3 — 13) 
(3 — 14) 
(3 — 15) 
(3 — 16) 
(3-17) 


Eq.(3-1) can be rewritten as 

(6 — 04) + K,(@ — 64) + K,(@ — 04) =0 | Ee 18) 
By substituting Eq.(3-14) into Eq.(3-18), we find the contro! angle (0) 

5 =—B™'[ A(6q— K,(0 — 64) — K,(0 — @,)) -F] (3 — 19) 


The control angle (6) is a function of the pitch angle (@) , desired pitch angle 

(8,), angular velocity (6), desired angular velocity (6,), desired angular acceler- 

ation (8,), the position feedback gain Kp , the velocity feedback gain K,, and the 

| matrices i.e., A, B, F. Since the thrust force T is included in F and B, the mag- 

nitude of the thrust force will thus affect the control angle 6 . K, and K, are ad- 
justable to obtain the desired response of the pitch angle @. 


IV. COMPUTER SIMULATION 


A. SIMULATION OBJECTIVES 

In literature, the application of forces has been limited to gravity, aerodyna- 
mic forces and thrust [Ref. 2]. No damping has been implemented. The thrust 
and aerodynamic forces can be found fairly accurately with standard test and 
design procedures such as Static firings to obtain thrust versus time profiles for 
the engine and wind tunnel measurements to determine the aerodynamic forces. 
Gravitational forces can be calculated from the knowledge of the’ missile’s posi- 
tion relative to the earth. The mass including fuel can be estimated from know- 
ledge of the missile’s weight before and after burnout (from the measurements), 
and by using a mathematical relationship (often linear) for the decrease in missile 
mass over the engine burntime. , 

In this work, the missile’s weight is assumed constant and the aerodynamic 
forces are zero. The deformations resulting from the structural flexibility have 
been assumed small and small control angle assumptions are used in the rigid- 
body controller design. 

The primary purpose of this study is to complete the simulation of the flexible 
missile in 2-D motion, where the bending effect is considered to be the only flex- 
ibility source and a rigid-body controller as developed in the last section is in- 
cluded. 


B. SOLUTION TECHNIQUE 

Many numerical integration techniques can be applied in the solution of the 
equations of the motion. The main consideration in the selection of the inte- 
gration technique is the size of the time step necessary to integrate the equations 
of motions with numerical accuracy and stability. 

The type of the equations of motion (2-20, 2-21) in this research permits the 
application of the Sequential Integration Method [Ref. 6]. The linear equations 
of small motions are integrated implicitly and the large motion solutions are ob- 


tained using explicit integration method. The implicit methods are effective for 
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linear systems with high frequencies and the explicit methods are effective for 
solving nonlinear systems with low frequencies. The implicit method is especially 
effective for linear systems having a wide range of frequencies of which only the 
lower frequencies are excited. The size of the time step need only be chosen to 
make the solution of the excited modes sufficiently’ accurate. Explicit methods 
are effective for large scale systems with low frequencies. Because of the low fre- 
quencies, the size of the time step that we choose for stability need not be small. 
In addition, the explicit method does not need iterative procedures, which are 


time consuming for non-linear systems. 


C. THE COMPUTER SIMULATION CODE 

A high level computer language, i.e.. MATLAB, was chosen to simulate the 
missile system. The MATLAB was designed for matrix operations. The simu- 
lation code was developed with modular MATLAB routines. 

The simulation code can be divided into three levels : LEVEL | (an overview) 
Separates the code into primary portions of INITIALIZATION, PLANT DE- 
SCRIPTION, INTEGRATION, SYSTEM CONTROL and OUTPUT. LEVEL 
2 facilitates several subroutines for LEVEL 1. A listing of MATLAB routines 
required for manipulation in LEVEL 2 is placed in LEVEL 3. 


D. SIMULATION RESULTS 
The computer simulation was performed with variable parameters which de- 
termine missile speed and controller bandwidth. These parameters include force 
iesand K,. 
The simulation outputs will be presented in large motions, small motions and 
control angle. The large motions are X,, Y, and @ and the small motions are v(0) 
and (0). The initial condition for all runs was the pitch angle of 45°. The sim- 
| ulation work was divided into two areas. First a simulation was performed for 
the rigid missile using the rigid-body controller. These results were used as a 
~ baseline for comparison with the results of the flexible missile using the rigid-body 
controller system. Second a simulation was performed for the flexible missile and 


rigid-body controller system. Only the trajectory control will be discussed in the 
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analyze of the dynamics of the flexible missile. The desired trajectory was speci- 
fied as 


45° when 0.1 <t<0.l 
@6=| 45° —112.5t when 0.1 <t<0.4 
0 when t>0.4 


The control angle is assumed to be limited to +/-10 degrees (small angle). 
This restriction puts a Saturation line to the control input. 

Fig.5 represents the graphical results of the desired and actual trajectory for 
the rigid missile and rigid-body controller. The force (T) and controller band- 
width (w,) were 30000 N and 3 rad/s respectively. The dark black line and 
dashed line represent the desired and actual trajectory, respectively. The control 
angle (0) is shown in Fig.6 and Figures 7-8 present the base positions (Xo Yo) in 
large motion. 

After presenting results for the rigid missile with rigid-body controller, the 
dynamic behavior of the flexible missile with rigid-body controller will be studied 
next. The force (T) and controller bandwidth (w,) again were 30000 N and 3 
rad/s, respectively. Fig.9 shows the desired and actual trajectory. The control 
angle is presented in Fig.10. After 0.1 sec, the missile initially needs a control 
angle which is less than 10 degrees. The effects of small motion can be seen on 
control angle clearly. The base deflection and slope are shown in Figures I1-12. 
There is no small motion effects on flexible missile between 0-0.1 sec because of 
zero control angle and no force components. After 0.1 sec, the small motions are 
excited and have an amplitude of 10-3. Figures 13-14 present the base positions 
of the flexible missile in large motion. The large motion behaves like the rigid- 
body motion, which implies that the coupling effect between the large motion and 
small motion is small. The dynamic behavior of the flexible missile is dominant 
by the large motion in this case since the bandwidth of the controller is low and 
Small motion is not significantly excited. 
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Figure 5. The desired and actual trajectory for the rigid missile with rigid-body 
controller (T = 30000 N, w, = 3 rad/s) 
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Figure 6. The control angle for the rigid missile with rigid-body controller (T = 
30000 N, w, = 3 rad/s) 
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Figure 7. The large motion position X for the rigid missile with rigid-body controller 


(T = 30000 N, w, =3 rad/s) 
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Figure 8. The large motion position Y for the rigid missile with rigid-body controller 
(T = 30000 N, w, = 3 rad/s) 
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Figure 9. The desired and actual trajectory for the flexible missile with rigid-body 


controller (T = 30000 N, w, = 3 rad/s) 
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The control angle for the flexible missile with rigid-body controller (T 


—soOv00 N,@e = 3 rad/s) 


Figure 10. 
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Figure 11. The small motion position v(0) for the flexible missile with rigid-body 
controller (T = 30000 N, w, = 3 rad/s) 
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Figure 12. The small motion position phi(9) for the flexible missile with rigid-body 


controller (IT = 30000 N, w, = 3 rad/s) 
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Figure 13. The large motion position X for the flexible missile with rigid-body con- 
troller (T = 30000 N, w, = 3 rad/s) 
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Figure 14. The large motion position Y for the flexible missile with rigid-body con- 
troller (T = 30000 N, w, = 3 rad/s) 
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Two kinds of tests were performed to further the study of the dynamic be- 
havior of the flexible missile. 

The first test was to study dynamics due to the increase of the control system 
bandwidth (w,). The desired and actual trajectory of the flexible missile with 
rigid-body controller without saturation on the control angle is presented in 
Fig.15. The force (T) is still 30000 N while the controller bandwidth (w,) is in- 
creased to 300 rad/s. As the controller bandwidth is increased, the gap between 
controller bandwidth and system natural frequency is smaller. The system will 
be unstable at controller bandwidths close to the system fundamental frequency 
which can be as high as 270 rad/s in this case. Note that the fundamental fre- 
quency of the simple beam with free-free boundary conditions was calculated as 
@,, =209.2053 rad/s, and the missile’s fundamental frequency will be increased 
due to the coupling between large and small motion. The control dynamics in- 
terfere with the structural dynamics. The higher control frequency with the 
bandwidth of 300 rad/s significantly excites the small motion of the structure. 
The graphs of control angle and the small motions and large motions of the mis- 
Sile base (Figures 16-17-18-19-20) show the unstable state. 

The bandwidth of the controller must be set much lower than the funda- 
mental frequency of the missile in order to use the rigid-body control. This implies 
that the pitch-angle response speed, which is determined by the control band- 
width is limited. To keep the response speed, the missile structure must be de- 
signed sufficiently rigid to possess a high fundamental frequency. On the other 
hand, the missile payload will affect the natural frequency of the missile structure 
with the heavier the payload, the lower the fundamental frequency of the missile. 
Therefore, the payload must be limited to achieve high-speed control response. 
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Figure 15. The desired and actual trajectory for the flexible missile with rigid-body 
controller without saturation on the control angle (YT = 30000 N, w, 
. = 300 rad/s) 
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Figure 16.  Tlte control angle for the flexible missile with rigid-body controller 
without saturation on the control angle (IT = 30000 N, w, = 300 
rad/s) 
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Figure 17. The small motion position v(0) for the flexible missile rigid-body con- 
troller without saturation on the control angle (T = 30000 N, w, = 
300 rad/s) 
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Figure 18. The small motion position phi(0) for the flexible missile rigid-body con- 
troller without saturation on the control angle (f = 30000 N, w, = 
300 rad/s) : 
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Figure 19. The large motion position X for the flexible missile rigid-body controller 
without saturation on the control angle (T = 30000 N, w, = 300 
rad/s) 
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Figure 20. The large motion position Y for the flexible missile rigid-body controller 
without saturation on the control angle (T = 30000 N, w, = 300 
rad/s) 
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Fig.21 presents the desired and actual trajectory of the flexible missile using 
rigid-body controller with a saturation (+/-10 Degrees) on the control angle. The 
force (T) and controller bandwidth (@,) are 30000 N and 200 rad/s, respectively. 
A better trajectory tracking is obtained, however the actual trajectory is oscillat- 
ing about the desired trajectory which is mainly due to the switching of the con- 
trol angle between the saturation lines. Fig.22 shows the control angle, and 
Figures 23-24 show the small motion displacement and slope of the base of the 
flexible missile. It was observed that the switching of the control angle between 
+/-10 degrees causes increased excitation of small motions. Figures 25-26 pre- 
sents the position of missile base in large motion where the large motion 1s no 
longer dominated by the rigid-body motion. 

The second test was to study the dynamics of the flexible missile due to an 
increase of the missile speed. To increase the speed of the missile, the thrust force 
was increased. Fig.27 presents the actual and desired trajectory of flexible mis- 
sile. The force was 60000 N and controller bandwidth was 3 rad/s. As seen from 
Fig.27, the pitch-angle response is determined directly from the bandwidth of the 
controller. The increased missile speed does not change the pattern of the pitch- 
angle response. The control angle of the flexible missile is shown in Fig.28. The 
control angle is affected by the increased speed. When the speed 1s increased, the 
control angle gets smaller because a smaller control angle generates a sufficient 
correction to the pitch angle. The increased missile speed has little effect on small 


motions (Figures 29-30). Figures 31-32 show the large motions. 
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Figure 21. The desired and actual trajectory for the flexible missile using rigid-body 


controller with saturation on the control angle (T = 30000 N, w, = 


200 rad/s) 
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Tigure 22. The control angle for the flexible missile using rigid-body controller with 


saturation on the control angle (T = 30000 N, w, = 200 rad/s) 
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Figure 23. The small motion position v(0) for the flexible missile using rigid-body 
controller with saturation on the control angle (T = 30000 N, w, = 
200 rad/s) 


39 


MISSILE CONTROL SYSTEM 
TIME(S) / 





0.01 


WY 
S 
S 


rr 
= 
a 


0.005 


© 


(ADNOILISOd HdOTS NOLLOW TIVIWS 


-0.015 


Figure 24. The small motion position phi(0) for the flexible missile using rigid-body 
controller with saturation on the control angle (J = 30000 N, ow, = 
200 rad/s) 
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Figure 25. The large motion position X for the flexible missile using rigid-body 
controller with saturation on the control angle (IT = 30000 N, @, = 
200 rad/s) 
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Figure 26. The large motion position Y for the flexible missile using rigid-body 
controller with saturation on the control angle (T = 30000 N, ow, 
= 200 rad/s) 
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Figure 27, The desired and actual trajectory for the flexible missile using rigid-body 


controller with increased speed (T = 60000 N, w, = 3 rad/s) 
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Figure 28. The control angle for the flexible missile using rigid-body controller with 
increased speed (T = 60000 N, w, = 3 rad/s) 
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Figure 29, The small motion position v(0) for the Mexible missile using rigid-body 


controller with increased speed. (T = 60000 N, w, = 3 rad/s) 
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Figure 30. The small motion position phi(0) for the flexible missile using rigid-body 
controller with increased speed (T = 60000 N, w, = 3 rad/s) 
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Tigure 31. The large motion position X for the flexible missile using rigid-body 
controller with increased speed (T = 60000 N, w, = 3 rad/s) 
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Figure 32. The large motion position Y for the flexible missile using rigid-body 


controller with increased speed (T = 60000 N, w, = 3 rad/s) 
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VY. SUMMARY 


A. CONCLUSIONS 

The development and computer simulation is presented for a bending flexible 
mussile with a rigid-body controller system moving in a 2-D coordinate frame. In 
this research, the dynamic model has been developed through the Equivalent 
Rigid Link System utilizing Lagrangian dynamics to obtain a type of system 
equations of motion suited for Sequential Integration Method that integrates 
large motion explicitly and small motion implicitly. The spatial finite element 
discretization of missile structure and the application of truncated natural modal 
responses provide an approximate solution. 

The analysis and simulation were performed to understand the dynamic be- 
havior of a flexible missile using a rigid-body controller. It was found that the 
controller bandwidth must be much lower than the fundamental frequency of the 
missile in order to use the rigid-body controller. The payload will affect the na- 
tural frequency of the missile structure i.e, when the payload is increased, the. 
system fundamental frequency will be decreased. The payload must then be lim- 
ited to achieve high-speed response. In order to increase the payload and main- 
tain high-speed control response, a flexible-body controller is needed. 


- 


B. RECOMMENDATIONS 
Areas which remain to be investigated include : 
|.Add the payload and aerodynamic effects to the model. 
2.Design and study the dynamic behavior of a flexible missile with a 
flexible-body controller in 2-D motion. 
3.Build a flexible missile in a laboratory scaie and obtain experimental data. 


4,Design and simulate a control system for flexible missiles in 3-D motion. 
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APPENDIX A. DERIVATION OF THE MODE SHAPE RESPONSE 
MATRIX FOR THE FLEXIBLE MISSILE 


In this study, the natural mode shape functions of a beam are used to repre- 
sent the flexural motion of the flexible missile. Only the first two mode shapes are 
used. The flexible missile is modeled as a continuous Euler-Bernoulli frec-free 
beam, neglecting shear deformation and rotary inertia effects. 


The bar (Fig.33) has system parameters: 


f(x,t) 
m(x), E(x) 


ea 


ax -— 


a 


a ae 
nel 





Figure 33. The Bar in Flexure 


y(x,t) = Transverse displacement at any point x and time t 
{(x,t) = Transverse force per unit length 

m(x) = The mass per unit length 

El(x) = Flexural rigidity 

E = Young’s modulus of elasticity 

I(x) =The cross-sectional area moment of inertia 

Q(x,t) = Shearing force 

M(x,t) = Bending moment 
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Fig.34 shows the free body diagram corresponding to a bar element of Icngth 
dx. | 








f(x,t) dx 
dQ(x,t) 
M(x. t) Q(x,t) dx 
ax 
Q(x, ¢) fen, t)> ne dx 
ax 


y(x,t) 


Figure 34. Free Body Diagram Corresponding to a Bar Element 


The force balance of the free body is 


200) Py(x,t) 


[O(x,.02) + ee] — Cl) yee — m(x)dx ——>— ae (A-— 1) 


- The moment equation of motion about the axis normal to x and y, ignoring 


the inertia torque associated with rotation, 


Sl 


[.M(x,t) + at dx | — M(x,t) + [O(x,t) + <0 axa ae 


+ flxdax & ay ia 3) 


Canceling appropriate terms, ignoring terms involving second order terms in 
dx and combining Eqs. (A-1) and (A-2), 


2 
a x j= M(x L289 (Ce 3) 
6t? 


Eq. (A-3) relates the bending moment M(x,t), transverse force f(x,t) and 
bending displacement y(x,t). 
The relation between the bending moment and the bending deformation is 


674(x,0) 


Peery — EI(x,f) 5 (4 — 4) 
Ox 


Inserting Eq. (A-3) into Eq. (A-4), we obtain the differential equation for the 


flexural vibration of a bar, 


2 
=F cen) ZH) + pes) = ates) as 


The bending moment and shearing force of the free ends (x=0, x=L) are 


Zero. 
ene yy OVE) 3 y(x,t) 
El(x) ——.— |.=o = 0 E(x) ——— x= = 9 (A= a6) 
Ox Cx 
Q a(x) Q "y(x,t) 
Gx IO) I s-0 = 9 Gy LEI be lame (A) 


Eqs. (A-6) and (A-7) are called natural boundary conditions. 


De 


Considering the free vibration characterized by f(x,t)=0, using separation of 


variables method and simplifying Eq. (A-5), 


F612) AO J = wm re9 (4 —8) 
ax dx? 


Simplifying the Eq. (A-8), for EI(x) =constant 





d’ 
<n (4-9) 
ae 

Wilere b> = a 


The boundary conditions require that, 
At x=0 (base) ; 











d’Y(-: 
ay =p ae 
dx 
aY(- 
¥(x) _ =0 (4 — 10) 
3 x=0 
ax 
At x=L (tip) 
a7 VC 
= =p = 0 
ax 
avy | 
Ce 


The general solution of Eq. (A-9), 
Y(x) = C, sin(Bx) + C, cos(Bx) + C3 sinh(Bx) + Cy cosh(Px) (A — 12) 


Taking the derivatives of Eq. (A-12), and substuting boundary conditions, 
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Y(x) = C,(sin Bx + sinh Bx) + C,( cos Bx + cosh Bx) (4 — 13) 
Solving Eq. (A-13) yields, 
Ges PL cosh BL = } (4 — 14) 


The first five consecutive roots of this equation are; 


B, = 0.0 


B, = 10.995574 
Eq. (A-13) is now written in the following form, 


Y (x) = C,( cos Bx + cosh Bx) + (sin Bx + sinh Bx) r=1,2 (4 — 15) 
where: 


sin B,L — sinh BL 


eae ee A A — 16 
" —cos B,L + cosh B,L ( 


The transverse displacement v(x) and slope M(x) can be represented in the 


following forms respectively, 


 Srare (4 —17) 
r=] 


(x) = (4 — 18) 
v(x) = a,(C,(cos fx + cosh B,x) +-(sin B,x + sinh Bx) 

+ a,(C>( cos fx + cosh Box) + (sin Box + sinh Bx) (4 — 19) 
(x) = a,(C,hyx + sinh B,x) + B,(cos B,x + cosh B,x)) 


ar a,(C>B>( — sin Bx + sinh B>x) - (cos Box + cosh B5x)) (4 a 20) 
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Substituting the boundary conditions into shape function gives, 


v(O) = 2a, Cy + 2a,C) 
D(0) = 2a), + 2a>f, 


The modal amplitudes a, and a,; 


l Copy C2 
1=(>c hee CoE )v(O) + GE?) 





= (0) + 0(0) 
y = av(0) + b&(0) 
a= F,Y,(x) + F3Y>(x) 


b = FrY\(x) + Fy Y2(x) 











Pe I Crp 
2C, " CE 
a5 

= Cig 

B 

————— 

; CE 

ls 

Fy= 

G 

E=2f,-2 a 


(A — 21) 


(4 — 22) 


(A ey 


(A — 24) 


(A — 25) 
(4 — 26) 


(A — 2am) 


(A — 28 


(A — 29) 


(4 — 30) 


(A 22 


(4 23a 


Substituting Eqs. (A-15), (A-16), (A-26), (A-27), (A-28), (A-29), (A-30), 
(A-31), (A-32) into Eq. (A-25) yields an expression for the transverse displace- 
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ment of a flexible missile as a function of the missile base nodal displacements, 
v(O) and (0), 


»(x,t) = [Fy(C,(cos Bx + cosh B,x) + (sin Byx + sinh B,x)) 

+ F3(C,( cos Bx + cosh f,x) + (sin B,x + sinh f,x))]v(0) 

+ [Fo(C,( cos B,x + cosh B)x) + (sin B,x + sinh B,x)) 

+ F3(C,( cos Bx + cosh Bx) + (sin B,x + sinh Bx))]®(0) (A — 33) 


This expression is differentiated twice to obtain v”(x), which is necessary for 


the calculation of the potential energy due to deformation and theoretical strain, 
v"(x,0) = (F,(C,B%( — cos Bx + cosh B,) + Br — sin Bx + sinh f,x)) 

“= a C585 ( — cos f,x + cosh B,) + B5( — sin B,x + sinh f,x)))v(0) 

fe (F>(C, By ( — cos B,x + cosh f,) + Bail — sin Bx + sinh f,x)) 

+ F,(C>P3(— cos Bx + cosh 85) + B5( — sin Bx + sinh B5x)))®(0) (a5) 


Now substitution of v(x) into the 3x1 deformation vector, d a yiglas the ox 2 


Shape function matrix, @, and the 2x1 nodal displacement vector, U, 


0 00 
aaov—| () |= 00 | (Ad — 36) 
(0) 
v(x) ab 


The shape function matrix 1s now in a form convenient for computer coding. 
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APPENDIX B. DERIVATION OF THE EQUATIONS OF MOTION FOR 
THE FLEXIBLE MISSILE 


The ERLS 1s used to separate the flexible missile’s motion into a small motion 
and a large motion. The large motion generalized coordinates are defined as Xp, 
Y,, theta (@)and the small motion generalized coordinates as U. The Lagrangian 
Dynamics are very helpful in the derivation of the equations of motion of com- 
plex systems. U represents a column vector which has two elements. First one is 
v(0) which is the nodal displacement of the flexible missile’s base and ®(0) which 
is the slope of the flexible missile’s base. We can represent the vector of the gen- 


eralized coordinates by Eq. (B-1), 
é = (eae U =[v(0) (0) J (B—1) 


These generalized coordinates will be used in the application of the Lagrange 
equations to the equations of motion for both small and large motion coordinate 
stems. The Lagrangian equations are defined in Eqs. (B-2) and (B-3). 

d-0KE , odKE , OPE 


at 0é, ] eé, at; oé, Fs; (i oro) ( ) 


2 ae) Ge, (B— 3) 
7 gU aU 
aU 


where €, is a component of € 
Using ERLS we can define the global position of the base position in terms 
of the local position vector (r), a deformation vector (d@), and a coordinate trans- 


formation matrix (W) with Eq. (B-4), 


R=W(r +4) (B—4) 


oY 


Continuing in the development of the equations of motions, we determine the 
derivative with respect to time for global positions in order to obtain the kinetic 


energy expressions. 


a 
¢ = 


R=W( +d)+ Wd (B — 5) 


The Kinetic and potential energy expresions follow. In this development we 
have seperated the large and small motion energy contributions and will present 


them seperately. 


[R™Rdm (B — 6) 


PE =~ JU"T CL Udx — [R"gdm (B-7) 


Redefining d in terms of the shape matrix @ (derivation of d is presented in 


Appendix A ), we can rewrite Eq. (B-6) as 


KE = > [\((We +4) + Wa} "LW +d) + Wd))am (B — 8) 
KE == [WTWr + 2 WTWOU + 2 WTWOU 
+U'¢' WwW Weu +2U'd WwW WweU + Ud WI WbU dm (B —9) 


Continuing in the development of the equations of motion, we first express 
the derivative of the Kinetic energy with respect to the time rate of change of the 
large motion coordinate €, Eq. (B-13), and then determine the time rate of change 


of this expression, Eq. (B-14). Since 


58 


aa Gy ee (B — 10) 
0 €3 
Wwe = = ae = Partial derivative w.r.t. X 
i < = = = Partial derivative w.r.t. Y 
ea ee = Partial derivative w.r.t. @ 
00 CC, 


W , = Partial derivative of W w.r.t. X, Y, @ 1=1,2,3 





, __——— iv, 
OKE =F I T = ie 4 pe — pe ow! WU 
0¢; 0¢; ie 


i ow? 


+ uTeT WT bU +2U7¢ WU ) dm (B11) 


i i 


Putting Eq. (B-12) into Eq. (B-11) and simplifying 








on = S = Wy, a =Wi (22 


7 ae = [7 wl +r WIWOU +7 WW OU 


i 
+7 WIWOU + U6 WI WOU + U6 TWTIWOU ) dm (B — 13) 
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al a.) =f wir +7 wii +7 WI WeU 
i 


+7 WIWOU + wrod +7 WW OU 

+7 WW 3oU + re +7 WW OU 

+7 Wivoutr wl 5 i 4 "wl WwoU 

+ Uo WIWOU + UTS WTWOU + Ub WT WOU 

ip etntwesi a ae yai x ey 

7 wai ) dm (B — 14) 


Finally we express the derivative of the kinetic energy with respect to the 


large motion position as ¢,. 


OKE 
ag, 


=(( wir +7 wlweu +77 WT Wo 
+7 wlwou +7 WW ou + Uo WIWobU 


+ Uo WiwoU + U6 WW ou + UTG™WTW OU ) dn (e— 15) 
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d ,@KE,_ @KE 


a -_? 7 [FTW +7 WT WOU + 2 WI WoOU 
i i 


+7 WW OU +7 Wi WOU + Uo WT WOU 


+ 2U'¢ WI WoU + U6 WIWOU) am (B — 16) 


Completing our development, the expressions for the potential energy are 
presented organized similarly to the previous material, 


Cie 


i |¢ + 60) Wigan = - fr’ Wigdm — luo’ wigam (B -17) 
i 


The small motion equations are handled in a similar way, and are easy to 
derive than large motion. 


ERE = |g WWE + 6WTWOU + OWT WOU)an (B— 18) 
eur 
O1ae {eee ae Fe Toe 
RE _ ig Ww Wr + 6 W WoU + 6 W WoU)dm (B— 19) 
aut 
d_( OKE ) _ (gly wr + 6 wir +o W WOU + o WTO 
4 (SKE ) = [gw Wr + oT WTI + OTWTWGT + OT WTIVG 

aut 

426 W WOU + 6 W WOU + 6° W WOU ) dm (B— 20) 
4_( OKE y_ ERE = f(g Wr + 6 WTWOU 

aU 


dt* 
aut 
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+26 WT WoU + 6" W WoU ) dm (B — 21) 
Potential energy, 


a = (0 CrUax — |e" wT Zam (B — 22) 


We can simplify these expressions further by the substution of the second 
time derivative of the coordinate transformation matrix W. The contributions to 
the coriolis and centripetal forces represented by IV. , termed the residual accel- 


erations. The contribution to the general forces is represented by IV, . 
; : 3 : 
W=W,+ >) We; (B — 23) 
i=] 
For large motion, 


j GTwt(s W €,+ W,)r + ews We + WOU + 2r WlIVOU 
Si 2] r r eI J2j r st 


+r'(LWe,+ W)W OU +r W, WOU + Ub Wid We, + WOU 


~~ 


+ 2U'b WIW6U + Uro WI WOU —F Wig =g"W bU)dm =F; (B24) 


For small motion, 


(oi wt wes Worse oTwis we4 WOU Tw wou 
(2, eit Wr + ow (2. pop t WOU + 2ro't bU 
oo j= 
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+o Ww WoU +T'Cru — ¢ Wg =F, (B — 25) 


Collecting terms and arranging the coefficients the two Lagrange equations 
can be written as the equations of motion for the large and small generalized co- 
ordinates. It is these equations which must be solved by computer simulation 
Code: 


For the large motion, 


fr? } yi W ir dm + 


3|{[r'wiwoUdm+ |. | fr’ wiwdam+ |- 
rs ae = , i ps lee 

jail JU'd' WW oUdn+ |! | [Ub Wi Wed 
fU'o Wi Wr dm 


F,,— [rT Wi Wr dm — fr" Wl W,pUdm — 2[r’ WI WeUdm 
—fu'o WiW,rdm —[U'¢'’ Wi W,oUdm -2[U'o' Wi WoUdm (B — 26) 
+ fr’ Wigdm + |g’ Wj oUdm 


For the small motion, 


: (lol WTW dm + [ow W dUdmyJé, + ([o7W Wodm]U + 
=| 


[loo W WodmU + [[b7 W.bam + [PCP dx]U = 


(ee oe rc 
lo’ W'gdm—|b' W' Wrdm + F, (B — 27) 


63 


By defining 

M,,{i/) = fr’ WiWw jrdm + [ri wiw iP Udm + fu 'o'WiW iP Udm 
+fU'¢ WiWrdm (i=1,2,3) (= 1,2,3) (B — 28) 
|r’ wl Wwodm + [U6 Wl Wodm 


My, =| |r Wi Wea + [U'S' Wi Woam (B — 29) 
[r' W5Wodm + [U'o' W3Wodm 


= i, (B — 30) 
M,,=|¢' W' Wodm (B — 31) 
G, =|26'W' Wodm (B — 32) 
K, =|! W.odm + [0 CT ax (B — 33) 


—_ 


ni = Fe; — Je WWF dm — fr WT W,pUdm — 2)F Wi Wed 


—fuTo wi Wrdm —fU'b WTW,.6Udm -2[U6 Wl WU 


+frTwigan+|g™WUdm (i =1,2,3) (B — 34) 
F.=[¢' W'edm — |b W Wyram + F, (B — 35) 
r i 
Fo= [Fo For Fos | (B — 36) 


The final equations of motion are written as 


Moe é + Vr U = Fy 


—- os —" — 


EM U 4G enemy 


oT 


Mang 


ss 





APPENDIX C. DERIVATION OF THE GENERALIZED FORCES OF 
THE EQUATIONS OF THE MOTION 


To derive the generalized forces the principle of virtual work is used. We first 
write the transformation between global coordinates and local coordinates, Eqs. 
(C-1) and (C-2) respectively, 


l l 0 0 l 
Poy— | cos(?) sin(@) |x (C — 1) 
Y Yy sin(@) cos(@) |} y 


l l QO QO l 
x | =| —X_ cos(@) — Ypsin(@) cos(@) sin(@) || X (C — 2) 
y Xg sin(@) — Y,cos(@) —sin(@) cos(@) |} Y 


From Eqs. (C-1) and (C-2) we can find x and y. 
x = —X, cos(@) — Yo sin(@) + X cos(@) + Y sin(@) (C — 3) 
y = Xp sin(@) — Yo, sin(@) + X sin(@) + Y cos(@) (C — 4) 


Taking derivative of the Eq. (C-3) and (C-4) doing necessary substations, we 


can find, 
dx = 5X cos(@) + SY sin(@) (C — 5) 
dy = —6X sin(@) + 6Y cos(@) (C — 6) 


Using the virtual work principle, 
dW, = [ —Tcos(6 + D(0))v(0) ](d8 + dD(0)) 
+ [7 cos(o + ®(0))]éx + [T sin(d + O(0)) |(dy + dv(0)) (C — 7) 


Based on small angle assumptions, 
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cos(5 + @(0)) = 1 (C 
and 
sin(d + (0)) = 5 + (0) (C =a) 
Doing the necessary operations and substitutions, 4 
511%, = (T cos(0)—TS sin(@)—T(0) sin(6))6X + (T sin() + T cos(0) + 
T(0) cos(6))5 ¥ + ( — Tv(0))5(0) + (75 + TH(0))5v(0) 
+ (-Tv(0)50(0)) | (C-9) 


From Eq. (C-10), one can write large motion generalized forces (F,), i.e., Bae 
(C-11) and small motion generalized forces (F)), i-e., equation (C-12), 


T cos(@) — Té sin(@) — TP(0) sin(@) 
F, =| T sin(@) + 76 cos(@) +7(0) cos(0) (C — 11) 
— Tv(0) 


—=—> 


F 7” alg — (C12) 


7 —Tv(0) 
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APPENDIX D. COMPUTER SIMULATION CODE 


EEEEEEEESEEEEEEEEEEEEEEEEESE DEFINITIONS OF VARIABLES SEEEEEEEEEEEEEEEEEEES EEE 


Wo? OF AP oO 019 OF CP OF CP OP AP OF CP CP a OP CP CP OP OP GO OP CP oP OP OF OP CP EO OP OP OP OF OF CP OP CP GP cP oP oP? oP oo 


aQ-7 =The coefficients for the numerical integration 

area =Cross sectional area of the missile (m**2) 

betal =Modal noce 1 

beta2 =Modal node 2 

count =Constant that defines the element positions in vectors. 
dtheta =Desired pitch angle(dgrs) 

dthetadot=Desired angular velocity (rd/s) 

dtheta2dot=Desired angular acceleration (rd/s**2) 

delta =Control angle (dgrs) 


E =Young's modules (pascal) 

epsilon =Parameter for integration coefficient 

fn =Small motion force coefficient matrix 
flexible=Constant that does the flexible part on and off 
fnn =Reformed small motion force coefficient matrix 
fq ‘-=Large motion force vector 


ftime =Finish time of the integration 
force =External force (Newton) 


gn =Gyroscopic coefficient matrix 

grav =Gravitational constant (m/s**2) 

gravvec =Gravity matrix 

h’ =Time step used in the integration 

izz =Moment of inertia about the z-axis (m**4) a 
jota =Parameter for integration coefficients 


kn =Stifness matrix 
k =Nunmber of the column of the matrix 
kp,kv =The position and velocity feedback gains 
] =Number of the row of the matrix 
lforce =Large motion force vector 
=The length of the missile (meter) 
imn =Constant (1m*n) 
mn =Reformed small motion inertia and coupling coefficient matrix 
mqn =Coupling inertia coefficient matrix 
mag =Large motion inertia coefficient matrix 
ma =Constant (mu*area) 
mpohi =Shape function matrix 
mphi2 =Second time derivative the shape function 


mu =Mass density of the missile (kg/m**3) 
n =Number of the integration 
phi =Small motion slope position 


phidot =Gnall motion slope velocity 

phi2dot =Small motion slope acceleration 

rigid =Constant that does the controller on and off 
=Local position vector 
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ge oP oP dP OP OP OP OP dP dP oP OP OP OP GP cP OP OP OP GO OP oP oP cP OP OP oP UP oP oP oP oP cP oP oP oP oP ce oP OP oP oP oP cP 


q =Large motion generalized position vector 
qdot =Large motion generalized velocity vector 
q2dot =Large motion generalized acceleration vector 
stif =Stifness matrix 


sforce =Small motion force vector 

slope =Trajectory path slope 

tempfirst=Temporary matrix for the Simpson's rule 
templast=Temporary matrix for the Simpson's rule 
tenmpl... 

temp6 =Temporary matrices for integration 

tenp?7... 

tenplO  =Temporary matrices for large and small motion coefficients 
theta =Large motion angular position (pitch angle-dgrs) 
thetadot=Large motion angular velocity (rd/s) 
theta2dot=Large motion angular acceleration (rd/s**2) 


time =Sec 

u =Small motion generalized position vector 
udot =Small motion generalized velocity vector 
u2dot =Small motion generalized acceleration vector 
um =Small motion position vector 

umdot =Small motion velocity vector 

Vv =Small motion deflection position (m) 

vdot =Sinall motion deflection velocity (m/s) 


v2dot =Snall motion deflection acceleration (m/s**2) 

wresid =Residual acceleration matrix 

wdot =Time derivative of the transfommation matrix 

wksi =The derivative of the transformation matrix with respect to the 
large motion generalized coordinates 


X =Large motion x-direction position (m) 

x =Local x coordinate 

xdot =Large motion x-direction velocity (m/s) v 
x2dot =Large motion x-direction acceleration (m/s**2) 
y =Large motion y-direction position (m) 

ydot =Large motion y-direction velocity (m/s) 

y2dot =Large motion y-direction acceleration (m/s**2) 
mnni,kni, 

kn2, fnl, 

fn2, £n33, 

gni,mgqgl, 


fql,fq2 =Temporary matrices for the calculation of the coefficients 
CLEXCZ, Cl; 

f2,£3,f£4=Shape function constants 

m22, £22, 

h2,b2,A, 

B,F,templ14, 

temp20 =Temporary matrices for the controller 

ml1,m12, 

m21,m22, 


FP HM A AP oF WP A cP OP GP OP AP HP OP GO OM cH OM OF AI® OF 1? CO CIO CH CH OP CH CP CH OP of OP OM cP OP OF OP OP CY OP OP OP cP OO oO CP oO OO CO 


f11,hl,bi=Constants for the controller 
=Time plot vector 


trolot 
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oP oP oP CP CP oP oP OP OP CP oP OP dP OP dP oP OP OP OP oP oP cP OP dP oP oP UP oP oP oP oP oP oP OP oP oP oP oP oP oP oP oP OP OP oP oP OP oP oP of 


% xplot =Large motion x-direction position plot vector % 
% yplot =Large inotion y-direction position plot vector % 
%$ thetaplot=Large motion angular position plot vector % 
% xdotplot =Large motion x-direction velocity plot vector % 
% ydotplot =Large motion y-direction velocity plot vector % 
% thetadotplot=Large motion angular velocity plot vector % 
$ vplot =Smnall motion deflection position plot vector % 
% phiplot =Small motion slope position plot vector % 
% vdotplot =Snall motion deflection velocity plot vector % 
% phidotplot=Small motion slope velocity plot vector % 
% Deltaplot =Control angle plot vector % 
% dthetaplot=Desired trajectory angular position plot vector % 
% dthetadeg, deltadeg, % 
% thetadeg =Angles converted from radian to degrees % 
$ % 
$S% % 


EEFEETEEEETEEEEEEEETEEEEEEEEEEEEEEEEEEEEEEFETEEEEBEEEEETEEEEEEEEEBE EES ESE 


ESEESTEEEEEEEEEEEEEEEEEEEEESEEE MAIN PROGRAM ESTEEEEETESEEEEEEEESEEEEEEEEEEE EES 
% LEVEL 1 % 
% This is the controlling program for flexible missile. It defines the % 
% variables, determines the coefficients for explicit integration and mode % 
% shape function, calls large and small motion coefficient matrices routines, % 
% large and small motion integration routines, controller routines, output $% 
% routines % 
% $ 
SEEEEEEEEEEEESEEEEEEEEELEEEEEEEEEEEEEEEESEEEEEEEEEEGEEESEEEEEEEEEEEEEEEEEEEEES 


{X,y, theta, xdot, ydot, thetadot, 1m, betal, beta2,E,1zz, force, . 

phi, delta, v, iota, epsilon, ftime,h, grav, ma, vdot, phidot,n, mn, flexible, .. 
count, mphi,mphi2, stif, gravvec, rlocal, x2dot, y2dot, theta2dot, v2dot,... 
phi2dot } = constants]; 
[n22, £22,h2,b2,A,F,B, dtheta2dot, dthetadot, dtheta,kv,... 

kp, rigid, tenp14, slope, temp20)=constantsila; 


[g, qdot, q2dot, temp6, wresid, wdot, wksi,w,... 
lforce, sforce, um, umdot, fnn, gn, kn,mn, countl,ml1,m12,m21,f11,hl,... 
b1)=constants2a (X, y, theta, xdot, ydot, thetadot, x2dot, y2dot, theta2dot) ; 


[u, udot, u2dot, mnni, kn1, kn2, f£n1, £n2, £n33,mnn, gnl,mqq,mqn,... 


fq,mqql, temp1, temp2, tenp3, tenp4, temp5, temp7,mgnl, tenps, fql, temp9, ... 
templ0, fq2) = constants2b(v, phi, vdot, phidot, v2dot, phi2dot) ; 


[tplot, xplot, yplot,... 

thetaplot, xdotplot, ydotplot, thétadotplot, vplot, phiplot, vdotplot, .. 

phidotplot, thetadeg, deltaplot, dthetaplot, dthetadeg, deltadeg] =... 
constants3 (f£time, h) ; 


{a0,al,a2,a3,a4,a5,a6,a7) = ccoef (epsilon, iota, h); 
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(cl,c2, £1, £2, £3, £4)=mphicons (1m, betal, beta2) ; 


for time = 0.0:h:ftime, 
(w, wot, wksi, wresid, lforce, sforce, um, umdot) = wmatderivvariant (o 
theta, X, y, xdot, ydot, thetadot, force, phi,delta,v, vdot, phidot, ... 
wresid, wdot, wksi,w, lforce, sforce, um, umdot) ; 


[mqq,mgn, fq, fq2)=lrgcof (wksi, rlocal,mphi, um, w, lforce, wresid,... 
wdot , undot, gravvec, ma, Im, n,mqq, magn, fq,mqql,temp7,mqnl,... 
tenp8, fql, temp9, temp10, 1mn, cl, c2, f1, f2,f3, £4, betal, beta2, fq2) ; 


1f flexible—1, 


(mn, gn, kn, £nnj)=smLcof (mphi, w,mqn,magq, fq, wdot, wresid,mphi2,stif,... 
gravvec, sforce, ma, lm,n, lm, fnn, gn, kn,m,mnnl, kn1, kn2, fnl, fn2, fn33,... 
mnn, rlocal,cl,c2, f1,f2, £3, £4,betal, beta2) ; 


(u, udot, u2dot ]=intsml (a0, al, a2, a3,a4,a5,a6,a7,u, udot, u2dot, ... 
tenpl, temp2, temp3, temp4, temp5, mn, gn, fnn, kn) ; 


end 


(q,qdot,q2dot] = intlrg(h, a0, a3, a6, a7,mqn,mqq, fq,q, qdot, q2dot,... 
temp6, u2dot) ; : 


(X, y, theta, xdot, ydot, thetadot, v, phi, vdot, phidot, x2dot,... 
y2dot, theta2dot, v2dot, phi2dot ]=chvar (qg, qdot, q2dot, u, udot, u2dot) ; 


if rigid—1, 
[delta, dtheta, temp20)=rigidcontrol (maq, fq2, force, theta,ml1,ml2,m21,m22,... 
fll, f£22,hl1,h2,bl1,b2,A,F,B, dtheta2dot, thetadot, dthetadot,... 
dtheta, kv, kp, tenp14, time, slope, temp20) ; 

end 


[tplot, xplot, yplot, thetaplot, xdotplot, count, deltaplot, dthetaplot] =gplot(... 
count, time, X, y, theta, xdot, tplot, xplot, yplot, thetaplot, xdotplot, thetadeg,... 
delta, deltaplot, dtheta, dthetadeg, dthetaplot , deltadeg) ; 


[ydotplot, thetadotplot, vplot,, phiplot, vdotplot, phidotplot, count J=gplot2(... 
count, ydot, thet adot, v, phi, vdot, phidot, ydotplot, thetadotplot, vplot, phiplot, .. 


vdotplot, phidotplot) ; 

count1=count1+1 

1£ counti=0, . 
force=0; 

end 

if countli=0, 


tl=xX; 
t2=y; 
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t3=theta; 
t4=xdot; 
t5=ydot; 
t6=thetadot; 
t/v=v; 
t8=phi; 
t9=vdot; 
t10=phidot; 
tlil=x2dot; 
tl2=y2dot; 
t13=theta2dot; 
t14=v2dot; 
t15=phi2dot; 
t16=delta; 


clear X y theta xdot ydot thetadot v phi vdot phidot delta 

clear fnl fn2 f£n33 fnn fq fql gn gnl kn knl kn2 1force m mn 
Clear mnl mgn mgnl mg mgql q qoot q2dot sforce tenpl tenp2 
clear temp3 temp4 tenp5 temnp6 temp? temp8 tenp9 u udot u2dot 
clear um umdot w wdot wksi wresid countl tenpl0 

Clear x2dot y2dot theta2dot v2dot phi2dot 

clear mll ml2 m21 f1 hl bl templl templ2 templ3 tempfirst templast 


X=tl; 
y=t2; 
theta=t3; 

» xdot=t4; 
ydot=t5; 
thetadot=t 6; 
v=t7; 
phi=t8; 
vdot=t 9; 
phidot=t10; 
x2dot=t11; 
ly2dot=t12; 
theta2dot=t13; 
v2dot=t14; 
_phi2dot=t15; 
delta=t16; 


(q, qdot,, q2dot, temp6, wresid, wdot, wksi,w, . 
lforce, sforce, um, umdot, fnn, gn, kn,mmn, Sant; ml1l,m12,m21, f1,h1, 
b1}=constants2a (x, y, theta, xdot, ydot, thetadot, x2dot, y2dot, .. 


theta2dot) ; 


{u, udot, u2dot, mmnl, kn1, kn2, fn1, fn2, £n33,mnn,gnl,mqq,man,... 
fq,mqql, templ, temp2, temp3, temp4, tenp5, temp7,mqnl1, teips, fql, temp9, ... 
temp10) = constants2b (v, phi, vdot, phidot, v2dot, phi2dot) ; 
me , 

end 


AZ 


[plot1, plot2, plot3, plot4 plot5, plot6)=al) 
: 5 ; plotl (tplot, xplot, yplot,... 
thetaplot, xdotplot, ydotplot, thetadotplot, dthetaplot) 
plot ,plot8, plot9, plot10, plot 11)=allplot2 (t 
plot, vplot, phiplot,... 
vdotplot, phidotplot, deltaplot) ee 


EESEEEEEEEEEEEEEEEEEEEEEEETEEEEEEEEEEEEEEEEEEEEEEEEEEEETEEETETEEEEETELEETETEBE 


% LEVEL 2 % 
% This function detemnines constant values and initial values of the % 
% variables and matrices. % 
% % 


EELEEEEEEEEEEEELEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEETEEEEEEEEEEEEEEEEEELEEEEEEETE ES 


function [(X,y, theta, xdot, ydot, thetadot, lm, betal, beta2,E,1izz, force,... 
phi, delta, v, iota, epsilon, ftime,h, grav,ma, vdot, phidot,n, lmn, flexible,.. 
count, mphi, mphi2, stif, gravvec, rlocal, x2dot, y2dot, theta2dot, v2dot, ... 
phi2dot] = constantsl 


area=0.0113097; 
dnr=4 3 
betal=4.712398/1m; 
beta2=7.853981/1m; 
count=0; 

delta=0; x 
F=2.0e+11; 
epsilon=0.25; 
flexible=1; 
force=30000; 
frime=0.4; 
grav=-9.8066; 
h=0.001; 
jota=0.5; 
1zz=0.000010178; 
mu=7861.05; 

r=10; 

Imn=Im/n; 
ma=area*m; 
phi=0; 

phidot=0; 
theta=pi/4; 
thetadot=0; 

X=0; 

xdot=0; 

v=0; 

vdot=0; 

ydot=0; 

yO; 
nphi=zeros (3,2); 
mohi2=zeros (3,2); 
stif=zerbds (3); 
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stif (3, 3)=E*1izz; 
gravvec=zeros (3,1); 


gravvec (3) =grav; 

rlocal=zeros (3,1); 

x2dot=0; 

y2dot=0; 

theta2dot=0; 

v2dot=0; 

phi2dot=0; 

PELEGEREELEREELELEEEEEELEREE DELETE TER ERE ELE REELEE SEEDERS ED ERE RE EEE TEESE ESE SS 
% 
% This function determines constant values and initial values of the % 
% variables and matrices. ; 
% 


TLSELEEEEFEEESEEEESESEEEEEEEESEESSSESSESEEEEEBETTEESTETEBETEEBEEBEEEEEBEBEEEBE 


function (m22, £22,h2,b2,A,F,B, dtheta2dot , dthetadot, dtheta, kv, ... 
kp, rigid, temp14, slope, temp20)=constantsla 


m22=0; 

£22=0; 

h2=0; ~ 
b2=0; 

h=0; 

F=0; 

B-0; 

dtheta2dot=0; 
dthetadot=0; 
dtheta=0.7854;%45 Degree 
kv=6; 

kp=9? 

rigid=1;.- 

temp14=zeros (2) ; 
slope=-2.6180; 


tep20=0; 

SEFEETEEEEESEEEEEEEEEEEELEEEEEEEEEEEEELESEEESEELEEELEEEEEEEEEEEESEEEEEELEEEEEEELESSE 
% % 
% This function determines constant values and initial values of the % 
% variables and matrices. % 
% % 


EEEEEEELEBEEETELEEEEETEEETETEBELETEELEEEEEEELESEEEEEEETETEEEEEEEEEEEEESEEEEE EGY 


function [q, qdot, q2dot, temp6, wresid, wdot, wksi,w,... 
lforce, sforce, um, umdot, fnn, gn, kn,mn, countl,mli,mi2,m21,f11,hl,... 
b1)=constants2a (X, y, theta, xdot, ydot, thetadot, x2dot, y2dot, theta2dot) 
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q-=zeros (3,1); 
qdot=zeros (3,1); 
q2cdlot=zeros (3,1); 
terp6=zeros (3,1)? 
wresid=zeros (3) ; 
wcdot=zeros (3); 
wksi=zeros (3,9); 
w=zeros (3) ; 
lforce=zeros (3,1); 
sforce=zeros (2,1); 
um=zeros (2,1); 
umdot=zeros (2,1); 
frn=zeros (2,1); 
gn=zeros (2) ; 
kn=zeros (2) ; 
mn=zeros (2) ; 
count1=0; 

q(1) =X; 

g(2)=y; 

q(3)=theta; 

qdot (1)=xdot; 

qdot (2) =ydot; 

qdot (3)=thetadot ; 
q2dot (1) =x2dot; be 
q2dot (2) =y2dot; 
q2dot (3)=theta2dot; 
ml1=zeros (2); 
ml2=zeros (2,1); 
m21=zeros (1,2); 
f11=zeros (2,1); 
hl=zeros (2,1); © 
bl=zeros (2,1); 


ESEEEEEEEEEEELE EEE EEEEEEEEEEEEEEEELEEEEEEEEEEEBEEEEEEELEEEEEEEEEEEEEEEETEEEE ED 


% % 
% This function determines constant values and initial values of the % 
% variables and matrices. % 
$ % 


SEEEEEEEELEEEEEEEEETE EE EEEETEEEEEEEEEEEEEETEEEEEEEEEEEEEELEEEEEEEEEEEEETEEEEEE 


function [u, udot, u2dot,mnnl,kn1,kn2, fn1, fn2, £n33,mnn, gnl1,mqq,maqn,... 


fq,maql, termpl1, termp2, temp3, temp4, temp5, temp7, mgni, tenp8, fql, temp9,... 
terpl0, fq2) = constants2b (v, phi, vdot, phidot, v2dot, phi2dot) 


u=zeros (2, 1); 
udot=zeros (2,1); 
u2dot=zeros (2,1); 
mnnl=zeros (2) 3 
-knl=zeros (2); 
kn2=zeros (2) 3 


75 


imE=Zeros (2,1); 
fn2=zeros (2,1); 
fn33=zeros (2,1); 
mon=Zeros (2) ; 
gnl=zeros (2) ; 
mogzeros (3); 
myn=zeros (3,2); 
fq=zeros (3,1); 
myql=zeros (3) ; 
Cempl=zeros (2,1); 
temp2=zeros (2,1) ; 
temp3=zeros (2,1); 
Cep4=zeros (2,1): 
terpS=zeros (2,1); 
tenp7=zeros (3) ; 
mqnl=zeros (3,2); 
Cenp8=zeros (3) ; 
fql=zeros (3,1); 
fq2=zeros (3,1); 
temp9=zeros (3) ; 
tenp10=zeros (3) ; 
u(1)=v; 

u(2)=phi; 

udot (1) =vdot; = 
udot (2) =phidot; 
u2dot (1) =v2dot; 
u2zdot (2) =phi2dot; 


SESEEEEEEEBEEEEEEEEEEEEEEEEEEEEEEEEEEBEBEEEEEEEEEBETEEEEEEEBEEEEEEESEEEEEBEE EE 
t 


% 
% 


EIEEEEEBETEEEEETEEEEEEEEEEETEEEE EEE EEEEEEEEBEEEEEBEBBEEEEEBEBEBEEE EEE BEET TEBE 


% 
% 
% This function determines the initial values of the output matrices. 
% 

% 


function (tplot, xplot, yplot,... 
thetaplot, xdotplot, ydotplot, thetadotplot,vplot,phiplot, vdotplot,... 


phidotplot, thetadeg, deltaplot, dthetaplot, dthetadeg, deltadeg) =.. 
constants3 (ftime, h) 


tplot=zeros (1, (ft ime/h) +1) ; 
xplot=zeros (1, (ftime/h) +1) ; 
yplot=zeros (1, (ftime/h) +1) ; 
thetaplot=zeros (1, (ftime/h) +1) ; 
xdotplot=zeros (1, (ftime/h) +1); 
ydotplot=zeros (1, (ft ime/h) +1); 
thetadotplot=zeros (1, (ftime/h) +1) ; 
. vplot=zeros (1, (ftime/h) +1) ; 
phiplot=zeros (1, (ftime/h) +1); 
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vdotplot=zeros (1, (ft ime/h) +1) ; 
phidotplot=zeros (1, (ftime/h) +1); 
thetadeg=0; 
deltaplot=zeros (1, (ftime/h) +1); 
dthetaplot=zeros (1, (ft ime/h) +1) ; 
dthetadeg=0; 

deltadeg=0; 


FEIEEEEEELETLESEEEEEEESEEEEETELEEEEEEEEEEEEEEEEEEEEEETEEEEETEEEEEELEEEEEEEEEEEEEES 
% % 
% Formlating the coefficients used in the sequential integration method. % 
% % 
EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEELESEEEEEEEEEESEEEEEEEEEEEEEEEEEEEEEEEEEEEES 


function [a0,al,a2,a3,a4,a5,a6,a7) = ccoef (epsilon, iota, h) 
a0=1/ (epsilon* (h*h) ); 

al=iota/ (epsilon*h) ; 

a2=1/ (epsilon*h) ; 

a3=0.5/epsilon-1; 

a4=iota/epsilon-1; ~ 

a5=h/2* (a4-1) ; 

aG=h* (1-iota) ; 

a7=lota*h; i 


ESEEEEEEEEEEELEEEEETEEEEEEEEEETETEEEEEEEEEEEEEEFEEEEEEEEEEEEEEEEEEEEEEEEEEEEEES 


3 % 
% Formulating the coefficients used in the mode shape function matrix. — % 
% % 


EETEEEEEEEEEEEEEEEEEEEEETEEEEELEEEEESETEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEES 


function [cl,c2, £1, £2, £3, £4]=mphicons (ln, betal, beta2) 


cl=(sin (betal *1m) ~sinh (betai*1m) ) / (-cos (betal*1m) +cosh (betal*]m)); 
c2= (sin (beta2* 1m) -sinh (beta2*1m) ) / (-cos (beta2*1m) +cosh (beta2* 1m) ) ; 
=2*beta2-2* (c2/cl1) *betal; 

£1=1/ (2*c1)+(c2*betal) / ( (c1*2) *e) ; 

£2=~(c2/(cl*e) ); 

f3=-(betal/ (cl*e) ) ; 


£4=1/e; 
SEEEEEEEETETETETEEEEEEEEEEEEELETEEEEEEEEEEEEEEEEEEEEEEEFETELEFEEETETEEEBEEET EG 
% : % 
%$ This function calculates the matrices which changes with tine. % 
$ 


% 
SEEEEEEEEEEEELEEEEEEEEEEEEEEEEEETLEEEEEEELEEEEEEELETEEEEEEEEETETEEELEGETEEEEEEES 
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-% to X, Y, theta: WKSI 








function [w, wdot, wksi, wresid, 1 force, sforce, um, umdot)=... 
winatderivvariant (theta, X, y, xdot, ydot, thetadot, force, phi, delta,v,... 
vdot, phidot, wresid, wdot, wksi,w, lforce, sforce, um, umdot) 


SEEELEFETTELIGETTTILETEEEETETETELAESELETELTETETELESELETEEEEEELELEEETETETETEBEEE 
%  Conputation of the transfonnation matrix w % 
EEBTELEEEEETETETEELEEEEEEEEEEEEEEEEEETELELEEEEELEEETEEETELETELELELEEEEESETEE ES ES 
w(l,1)=l1; . 

w(2,1)=X; 

w (2,2)=cos (theta) ; 

w (2, 3)=-sin (theta) ; 

Ww (3, 1)=y; 

w (3,2) =sin (theta) ; 

w (3, 3) =cos (theta) ; 


TELTEBEELEEEEEEEEEEESEESSIFEEEEEESEESESEEESEIFEETEEEEESEISEESEEEEEESEEEEEESEESS 
% Computation of the time derivative of the transformation matrix WOOT % 
FEEEEEEEEEETELETEEEEEEEEEEEEETEEEESEESETEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE SBE 
wdot (2, 1)=xdot; 

wdot (2, 2)=—-thetadot *sin (theta) ; 

wdot (2, 3)=-thetadot *cos (theta) ; 

wdot (3, 1)=ydot; 

wdot (3, 2)=thetadot *cos (theta) ; 

wdot'(3, 3)=-thetadot*sin (theta) 7 


EEEELEEEEEEEEEEEETEEEETEEEEEEEEEEEEEEEEEEEEELEEEVEEEEEEEEEEEEEBEEEEEEEEEEEEE SS 

$ Computation of the derivative of the transformation matrix with respect % 
$ 

SEEEEEEEEEEEETEEEEEEETLEEEEEESEEEEETEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE 

wksi (2,1)=1; : 

wksi (3, 4) =1; 

wksi (2, 8)=-sin (theta) ; 

wksi (3, 8)=cos (theta) ; , 

wksi (2,9)=-cos (theta) ; 

wksi (3, 9)=-sin (theta) ; 


SEEEEEEEEEEEEEEEEEEEEEEEEEELELEEEEEEEEEEESEETEEESEEEELETETLELEEEEEEEETLETELEEEEES 
% .Compttation of the the residual acceleration matrix WRESID % 


FEEEEEEEEEELELE EEE EEETEEEEELEBEEEELEEEESSSEESSEEEEEEEEEEEEEEEEFEEEEEEEEESEEFESF 
wresid (2, 2) =~ (thetadot*2) *cos (theta) ; 

wresid(2, 3) =(thetadot*2) *sin (theta) ; 

wresid (3, 2) =~ (thetadot*2) *sin (theta) ; 

wresid (3, 3)=- (thetadot%*2) *cos (theta) ; 


EEEETEEELELELELEEEEEEEEEEEEEEEEEEEEEEEEETEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEESF 


% Computation of the large motion force matrix % 
SEEEELEEETELETLEEEEEELETEEELELETELEETEEEIEIAEEISISESSESSEEESESERSESSESRSSSSRSREF 


1 force (1) =force* (cos (theta) -delta*sin (theta) -phi*sin (theta) ); 
1 force (2) =force* (sin (theta) tdelta*cos (theta) +phi*cos (theta) ) : 
lforce (3)=-force*v; 
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FEEEEEEEEEEEEEEEEELEEELEEEEEELEEEEELEEEEEEESEEYESEESSEESEEEEEEEEESESESESEESEES 
% Computation of the small motion force matrix % 
EEETELETTELEEEEEEEEBEEEEEBEEETELEBEEEEEEEEESEEEEEEBESEEEEEESEESEEEEESEEEEEEEEES 
sforce (1) =force* (phi+delta) ; 

sforce (2) =v* force; 


TEEEETETEEEEELEBEEEEELEEEEEEELEEEEESEBEEEEEEESESESEEEEEEEEEEEEEEESEEESESEEEYSSF 
% Computation of the small motion position vector % 
SELEEBEEEEEEEELEEEEEEEEEEETEEEEBEEEEEEYESEEEEEEBEEEEEEEEEEEEEEEESESEEEEEEEESESF 
un (1) =v; 


um (2) =phi ; 
EEETEEETEEETEEEEEEEEESEEEEEELEETEEEEEEEEEEEESESEEEEEEEEESEEEESEEEESEESEESEEESS 
% Conputation of the small motion velocity vector % 


SETETEEEEEEEETEEEEELEEEEESEEEEEEELEEESESESESESESEEEEEEESEEESESESESEEEEEEEEEEEF 
umdot (1)=vdot; 
uindot (2) =phidot; 


SEEETEEELEEELEEEEEEEEELEEEEEEEEEEEESESEEEEEEEEEEEEEESEEEEEEEEEEEEESEEEEELEELE SS 
$ -—e 
% This function calculates values for the large motion coefficient matrices% 
% % 
SEEEEEEEEEEEEEEEEELEEETLEEEEEEEEEEEELEEETEEEEEEEEEEESEEEEELEEEEEEEEEEEEEEEEEEETS 


ed 


( 


function (mqq,man, fa, fq2)=lrgcof (wksi, rlocal,mphi,um,w, lforce,... 
wresid, wdot, umdot, gravvec,ma, lm,n,mqg,mgn, fg,mqgl, temp7,mgnl,... 
teip8, fql, temp9, temp10, lmn,cl,c2, f1, £2, £3, £4, betal, beta2, fq2) 


ESEEEEEEEEESEEEEEEEEEEEESEEEEEEEEEEEEEEETETEEEEEETEEEEEEEEEEEEETEEEEEEEEEEEEEES 


% Creating the coefficient matrix for MOO $ 
EEBESEEEEEEEEELESEEEELEBEEEEEEGEEEBESEBSEEEEEEEEBEEEEEEBEEEEEBESEBEEEESELSEEEES 
for i=1:3, 

templ=wksi (3, 3* (i-1) +1:3*4); 

for j=1:3, 


temp2=wksi (:, 3* (j-1)4+1:3*}4); 
temp7=tenpl '*temp2; 
maql (i, j) =quadmat (‘astringl11’,1,1,n,1mn,ma,cl,c2,betal,beta2, £1, £2,£3,... 
£4,mphi, 0,0, rlocal, temp1, tenp2, um, 0,0,0,0,temp7,0,0,... 
0, 0,0); 
end 
end 


ncgq=ma* (mqql) ; 


TECTUELE TELE TELEVTETTELETELEEEEEEEEEEVEEELEEEETELESELEEEBELESEBETEBELESEEEEESEE 


% Creating the coefficient matrix for MON % 
FEEEEEEEEEEEEEEESEBEEBEBEEEEEEESESEEEEEEEEEESESESESEESESESEBEEEESESEEEEESEEESS 


for i£1:3, ; 
temp=wksi (3, 3* (i-1) +1:3*1); 


le 


tenp8=terp' *w; 

mgnl (1, :)=quadmat (‘astring22',1,2,n, lm,ma,cl, c2,betal, beta2, £1, £2, £3, ... 
£4,mphi, w, 0, rlocal, 0, 0, um, temp, 0, 0, 0, 0, temp8, 0,0, 0,0); 

end 

mgnsma* (mgni)z 


LZEETEEBEEELEBEESEEEEEEEEBEBEEEEEESEEBESSEEEESEESEEEBEEBEESEE BETTE EBUEEEVE EVES 
% Creating the coefficient mayrix for FQ % 
SEFSEEEBEEEESEFEEEBESESEEBEEEEEEEBEESESEBSEEEEEEBEEEEEEEBEEEBUEBERETEEEEBEEEEGE 
tComputation of the FQ coefficient 
for i=1:3, 
temp=wksi (:, 3* (i- 1)+1: 3*i)3 
tenp9=teamnp' *wresid; 
tenp10—tenp '*wdot ; 
fql (1) =quadmat (' astring33',1,1,n, 
mohi,0,wresid, rlocal,0, 0, um, temp, 
tenp10, 0,0); 


Im/n,ma,cl,c2,betal,betaz2, fulpie2, £3,£4,..-.- 
wdot, umdot, gravvec, 0,0,tenp?, ..- 


end 
fq2=ma*fql; 
fq=fq2tlforce; 


FSEETEEEEESEELEEEEEEEEEEEEEEELEEEEETEEEEEEEEETEEEEEEEEEEBEEEEEEEEEEEEEBEEEBETEE? 
% 
% This function’calculates values for small motion coefficient matrices. 
% 
EVLEEESEEESELELEEEEEEEEEEEEEEEEEEEEEEEEEEESETEEEEEEEEESEFEEEEEBETEEEEEBESEEEEES 


ar gio 


oo 


. function {mn,gn, kn, fnn)j=smlicof (mphi, w,mgqn,mqgq, fg, wot, wresid, mphi2, ... 
stif, gravvec, sforce,ma, lm,n, Im, fnn,gn,kn,mn,mni,kn1,... 
kn2, £fnl, fn2, £n33,mnn, rlocal,cl,c2, £1, £2, £3, £4, betal, beta2) 


TESESEEEEEEEEESEEEEBESELEEELEEEEEEEEEEEEEEEEEEEEEEEEEEEEEGESELETEEEEEEEEEESESS 
% Creating the coefficient matrix MN % 
FESTEEETELEEEEEESEEESESEEEEESEEEEEEEEEESEEEESELEETESEEEEETEEEEEEEEEEEEEEEEEEET 
mnn=quadnat (‘astring44',2,2,n, lmm,ma,cl,c2,betal,beta2, £1, £2, £3, £4, ... 

mphi,w, 0, 0,0, 0,0, 0, 0, 0,0, 0,0, 0,0, 0,0); 
moni=ma*mnn; 
mn=mnn1—-mgn ' * inv (mqq) *mgn; 


FEETELELEEEEEEEEELELEEEEESESEESEFEEESESELEEEEEEEEEEEEEEEEEEEEFEEEEEEEEEEEEEEES 
% Creating the coefficient matrix GN % 
FEEBELEEEEEESEEESEEEGEEEEEEEEEEEEESESEEEESEEEEFEEEELEEEEEESESEEEEESEEEEESEEGEES 
gnl=quadmat (‘astring55',2,2,n, lm,ma,cl,c2,betal,beta2, f1,£2,£3,£4,... 

mphi,w, 0, 0,0, 0,0, 0, wdot, 0, 0,0, 0,0, 0,0, 0) ; 
gn=ma*gnl ; 


TEFTEEEEEEFTEEETESETEEEEEEEEEEEEEEGETETEESEEEEEEEETEEEEEEEEEEEEEEEEEEEEEEEEEEEET 
% Creating the coefficient matrix KN % 
FEEEEEETEGEFEEEEEEEEEFEEEEEEEEETEEEBEEEEEEEEEEEEEEGEEEEEBEBEBEEEEEEEEEEEEESGS 
kn1=quadmat (‘astring66', 2,2,n, lm,ma,cl,c2,betal, beta2,f1,f2,£3,f£4,... 
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mpohi,w,wresid,0,0,0,0,0,0,0,0,0,0,0,0,0, 0); 

kn2=quadmat ('astring77', 2,2,n, lm,ma,cl,c2,betal,beta2, £1, £2,£3,£4,... 
0,0,0,0,0, 0, 0,0, 0, 0,0, 0,0, 0,0, mphi2, stif) ; 

kn=ma*kni+kn2; 


EEELEEEEEEEEEEEEEEEEEEEEEEEEELEBEEEEEEEEEEEEEEEREEEEEEEEEEEEEEEEEEEEEEEREREE SE 

% Creating the coefficient matrix EN $ 

SELEEEESEEEEEEEEEEESESESEEESEBEEEEEEEEEEFESEFEEEEEEEEEEEEEFEEEEEEEEBEBBEEEEEESG ST 

fnl=quadmat ('astring88',2,1,n, Imm,ma,cl,c2,betal, beta2, £1, £2, £3, £4,mphi,w, ... 
wresid, rlocal,0,0,0,0,0,0,gravvec, 0,0,0,0,0,0); 

f£n33=ma* (fnl) + sforce; 

fnn=f£n33-mgn ' * inv (mag) * fq; 


ELEESEEEEEEEEEEEEBEESFEEEEEEEEBEEESESEEEEEEEEEETEEEEBEEEEEEEEEEVETEEEBETEBE ETE 


% 
% This function implicitly integrates the small motion equation. 


% 
SEEBEEEBESESEEEEESTSEEBEESEEEEEETEBEBBSEEBSVSEBEBEEBSEVEEBEBEBUEBVEEEBETEE VEE ETEEES 


49 


cP ol? 


function {u, udot, u2dot J=intsml (a0, al, a2, a3,a4,a5,a6,a/,u,udot,... 
u2dot, terpl1, terp2, temp3, tenp4, temp5, mn, gn, fnn, kn) 


terpl=a0*uta2tpdot+a3*u2dot ; 
temp2=al *uta4*udott+a5*u2dot; 
temp3=mn*tenpl ; 
tenp4=gn*temp2 ; 
tenpl=fnnttemp3+ttenp4; 
temp2=u; 
temp5=knta0*mntal*gn; 
templ=tempS\tenpl; 

u=templ ; 

terp3=u2dot; 

u2dot=a0* (u-temp2) -a2*udot-a3*tenp3; 
udot=udot+a6*temp3ta7*u2dot ; 


FEEEEEETELETETETEEEEESEEEEEEETELEEETELELELEEESEEEEETEEESEEEEEEEEEETEEEEEEE EES 
% 

% This function explicitly integrates the large motion equation. 

% 2 
EETEEEFELEEETESELEEEEEEETETESESEEEEEEEEEEEEEEEESEEEEEEEEEEEESESEEEEESEGEEEEEES 


a9 ¢.9 


ci 


ao 


function [q, qdot, q2dot ]=intlrg (h, a0, a3, a6, a7,mon, mqq, fq,q, qdot,... 
q2dot, temp6, u2dot) 


q=-gth*qdot+a3*q2dot/a0; 


qdot=qdot+a6b*q2dot ; 
temp6=temp6+mgn*u2dot ; 
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Eqn=£q-terp6; 
q2dot=inv (ingg) * £gn; 
qdot=qdot+q2dot *a7; 
q=qtq2dot/a0; 


TEEEEEEELEEELELEEEEEEEEEEEEEEEEFEEELEEEEEEEEEEEEEEEEEEELEEEEEEEEEGEEEEEEESSE os 
% 
% This function updates the time dependent values. 


% 
EEEEEEEEETEEETESETEELETEEETEEETEEETEEEEEEEEEEETEEEEEEEEETEEEEVEGEEEEEBEEEEBEE GGG 


a 


ci? ov? 


function [X, y, theta, xdot, ydot, thetadot, v, phi, vdot, phidot, x2dot,... 
y2dot, theta2dot, v2dot, phi2dot ]=chvar (q, qdot, q2dot, u, udot, u2dot) 


X=q (1)? 

y=q (2); 

theta=q (3); 
xdot=qdot (1); 
ydot=qdot (2); 
thetadot=qdot (3) ; 
x2dot=q2dot (1); 
y2dot=q2dot (2) a > 
theta2dot=q2dot (3); 
v=u (1); 

phi=u (2); 
vdot=udot (1) ; 
phidot=udot (2) ; 
v2dot=u2dot (1); 
phi2dot=u2dot (2) ; 


EFEEEEEELELEEEEELEFEEEEETEEEEEEEEEEEEEEEEEETEEEEEEEEEEBEEESEEEEETEEEEEEEESEET ES 
$ z 
% This function controls rigid or flexible missile with rigid-body s 
% controller. The control angle(delta) is limited to 10 degrees. ‘ 


% 
ELEEEEEEEELELEEETEEEEBEEELEEELELEEEEELELEEEEEEETEELEBEEELELEEEEEGELEEEVEGEEES 


9% 


o 
ow 
o 


function [delta, dtheta, tenp20)=rigidcontrol (mqq, fq2, force, theta,m11,m12,... 
m21,m22, £11, £22,h1,h2,b1,b2,A,F, B, dtheta2dot, thetadot, dthetadot,... 


dtheta, kv, kp, temp14, time, slope, temp20) 


if (time >= 0) & (time < 0.1), 
dtheta2dot=0; 
dthetadot=0; 
dtheta=pi/4; 
temp20=dtheta; 

elseif time >= 0.1, 
dtheta2dot=0; 
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dthetadot=s lope; 
dtheta=temp20+slope* (t ime-0.1); 
end 


m11(1,1)=mqq(1, 1)? 
mil (1, 2)=mqq(1, 2); 
ml 1 (2, 1)=mqq(2, 1)? 
m1 (2, 2)=mqq (2, 2); 
m2 (1)=mqq (1, 3)? 
ml 2 (2) =mqq (2, 3)? 
m22=mqq (3, 3)? 


£11 (1) =fq2 (1); 
£11 (2) =fq2 (2); 
£22=fq2 (3); 


hl (1)=force*cos (theta) ; 
hi (2) =force*sin (theta) ; 


bl (1)=-force*sin (theta) ; 
b1 (2)=force*cos (theta) ; 
tenp14=inv (m11) ; 
A=m22-m21*temp14*m12; 
F=f£22-m21*tenpl4* (£11+h1) ; 
B=m21*temp14*bl; 


delta=-inv (B) * (A* (dtheta2dot-kv* (thetadot-dthetadot) -kp* (theta-... 
dtheta) ) -F) ; 


if delta >= 0.174532 , %10 degree 
delta=0.174532; . 

end 

if delta <= -0.174532, 
delta=-0.174532; 

elseif (delta < 0.174532) | (delta > -0.174532), 


delta=delta; 
end 
TEELEEEEETEEEEEEETTEEEEETEETEEEEETEEEEEEEETEEETEEETEEEEETEEEEEEEEEETELESEEEEE ES 
% % 
% This function calculates the elements of matrices for output % 
: : % 


EEEELEEEEEEEEEEETEEETETETETELESESETEEESETETETESETEEEEEEESETELELEEELEEEEEEEEEEE 


function [tplot, xplot, yplot, thetaplot, xdotplot, count, deltaplot,... 
dthetaplot) =gplot (count, time, xX, y, theta, xdot, tplot, xplot, yplot,... 


s 
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thetaplot, xdotplot, thetadeg, delta, deltaplot, dtheta, dthetadeg, ... 
dthetaplot, deltadeg) 


count=counttl 
thetadeg=(180*theta) /pi; 
dthetadeg=(180*dtheta) /pi; 
dle ltadeg= (180*delta) /pi; 
tplot (count) =t ime; 

xplot (count) =X; 

yplot (count) =y; 

thetaplot (count) =thetadeg; 
xdotplot (count ) =xdot ; 
deltaplot (count) =deltadeg; 
dthetaplot (count) =dthetadeg; 


SEBEEEETEEEEETETEEEETEEETEEEELEFEEEETEEEEBEEEEBEEEBEEEEEEEEEEEBETETEEEEBETETESE 


% % 
% This function calculates the elements of matrices for output % 
% 


% 
EEEEEEEEEEEEEEETEEEEEEEEEEEEEEETEEEETEEEEEEEEEEEEEETEESEEEEEEEEEEEEETEEEEEEEEESES 


function {ydotplot, thetadotplot, vplot, phiplot, vdotplot, phidotplot, countJ=... 
gplot2 (count, ydot, thetadot, v, phi, vdot, phidot, ydotplot, thetadotplot,... 


vplot, phiplot, vdotplot, phidotplot) 


* ydotplot (count) =ydot ; 
thetadotplot (count) =thetadot; 
vplot (count) =v; 


phiplot (count) =phi; 

vdotplot (count) =vdot ; 

phidotplot (count) =phidot; 
FEETETEEEEEEEEETETEEEEEEEEEEEEEEEEEEEEEEEEEEEELEEEEEESESESEEEEEEEEEEEEEEEEEEES 
$ , % 
% This function plots the output graphs. % 
+" % 


FEEEEEEEEEEEEEEEEEEEEEEEELEEETELE EEE EBEEEEETEEEEESEEEBEEEESEEEEEEEEEEEEEBEEESBS 


function [plot1,plot2, plot3, plot4, plot5, plot6)=allplotl (tplot, xplot, yplot,... 
thetaplot, xdotplot, ydotplot, thetadotplot, dthetaplot) 


plotl=plot (tplot, xplot) , pause 

plot2=plot (tplot, yplot) , pause 

plot3=plot (tplot, thetaplot, '—-', tplot, dthetaplot, '+'),pause 
plot4=plot (tplot, xdotplot) , pause 

plot5=plot (tplot, ydotplot) , pause 

plot6=plot (tplot, thetadotplot) , pause 
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TEEEEESEEBEEEEEEEBEEBEEEEEEEEEEEEBEEBEEBEEBEEBEEEEBEBEEEEVEEELEVEEEEEEEEEEEE EEE 


% % 
% This function plots the output graphs. % 
% % 


EEEEEEEEEEEEEEEEEEEEEEEEEELEEEEEEETEEELETEEEBEBEEETEEEEEEEELELELEEEEEEETEEEEED 


function [plot7,plot8, plot9,plot10, plot11)=allplot2 (tplot, vplot, phiplot,... 
vdotplot, phidotplot, deltaplot) 


plot7=plot (tplot, vplot) ,pause 

plot 8=plot (tplot, phiplot) , pause 
plot9-plot (tplot, vdotplot) , pause 
plot10=plot (tplot, phidotplot) ,pause 
ploti1=plot (tplot, deltaplot) ,pause 


EELESSSEEEEEEEELELELELEEEEEEEETEEEELEEEEELELELELEEEEELEEEEEEEEEEEEEELEEEEEEEES 
% IEVEL 3 % 
% This function will integrate the coefficient matrices of small and large % 
% motion over the length of the missile using Simpson's integration method. % - 
% Z % 
EELESSLEEELELETEEEEEEELEETELEEEEEEEEEETELEEEEEEEEEEELEELELEEEEELELEEELELELELESS 


function aa = quadmat (F, 1, k, n, int, ma,cl,c2,betal,beta2,... 
f1,f2, £3, £4,mphi,w, wresid, rlocal,templ, teanp2, um, temp, wdot,umdot,... 
gravvec, temp7, temp8, temp9, temp10,mphi2, stif) 


templi=zeros (1,k*n); ‘ 
templ2=zeros (1, k); 

temp13=zeros (1,k); 

tempfirst=zeros (1,k);3 

templast=zeros (1,k); 

aa=zeros (1,k); 


x=0; 

for i=l:n, 
x=x+int; 
templl (:,k* (i-1) +1: i*k) =feval (F,x,cl,c2,betal, beta2,f1,f2,f3,f4,... 
mohi,w, wresid, rlocal,temp1, temp2, un, temp, rlocal,wdot, undot,... 
gravvec, temp7, temp8, temp9, temp10, mphi2, stif) ; 


end 


for m= 1:k, : 
for i = 2:2: (n-1), 
templ2(:,m) = templ2(:,m) + templl(:, k* (i-1) +m); 
end 
end 
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for p= 1:k, 

for j = 3:2:(n-l), 
temp13(:,p) = templ3(:,p) + templl(:, k*(j-1)4p); 
end ; 
ent 


tenpfirst (:,1:k)=templl (:,1:k); 
tenplast (:,1:k)=temp11 (:,n*k- (k-1) :n*k) ; 
aa=(int/3) * (tempfirst+4*temp13+2*tenpl2+ttemplast) ; 


oo 


SEVEEEEEEEEEETEEEEEEEEEEEEEEETIBESETEEEEELEEBEEEETEBEEEEEEEEEBEEEEEEEETEEEBEEES 


2 % 
% This function creates the product of matrices for integration of MQ % 
% coefficient matrix. % 
% $ 
SESEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEES 


function yll=astring]1 (x,cl,c2,betal, beta2, f1, £2, £3, £4,mphi,w,... 
wresid, rlocal, temp1, temp2, um, temp, rlocal, wdot, umdot, gravvec, temp7,... 
tarps, temp9, teinp10, mphi2, stif) 


mht (3, 1)=f1* (cl* (cos (betal*x) +cosh (betal*x) ) +sin (betal*x)+... 
sinh (betal*x) ) +£3* (c2* (cos (beta2*x) +cosh (beta2*x) )+... 
Sin (beta2*x) +sinh (beta2*x) ); 


| mphi (3, 2) =£2* (cl* (cos (betal*x) +cosh (betal*x) )+sin (betal*x)+... 
Sinh (betal*x) )+£4* (c2* (cos (beta2*x) +cosh (beta2*x))+... 
Sin (beta2*x) +sinh (beta2*x) ); 


rlocal (1)=1; 
rlocal (2) =x; 
4 
yll=rlocal' *temp7*rlocal+rlocal'*temp7*mphi*umtum!' *mohi'*temp7*... 
npohit*tumtum!’ *mphi '*temp7*rlocal; 


2° TESEEEEEEEEEEEEEEEEEEEEEBEEBEEEEEEESEBESSEEEEEEEEEEEEEEEELEBEEEBELEEEBEEEEEE 
2 % 
% This function creates the product of matrices for integration of MN % 
% coefficient matrix. % 
% % 
SITEESESEEEEEEEEEEEEEESEEEEEEESEEEEEBEEEEEEEBETTETEEEBETEVETEEEBEEELEEEEEEELETE 


function y22=ast ring22 (x,cl,c2,betal, beta2, f1, f2, £3, £4,mphi,w,... 
wresid, rlocal, terpl1, temp2, um, temp, rlocal, wdot, umdot, gravvec, temp7,... 
tenp8, temp9, temp10, mphi2, stif) 


ngohi (3, 1) =£1* (cl* (cos (betal*x) +cosh (betal*x) ) tsin(betal*x)+... 
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sinh (betal *x) ) +£3* (c2* (cos (beta2*x) +cosh (beta2*x))+... 
sin (beta2*x) +sinh (beta2*x) ); 


moh i (3,2) =£2* (c1* (cos (betal*x) +cosh (betal*x) )+sin (betal*x)+... 
Sinh (betal *x) ) +£4* (c2* (cos (beta2*x) +cosh (beta2*x))+... 
sin (beta2*x) +sinh (beta2*x) ); 


rlocal (1)=1; 
rlocal (2) =x; 


y22=rlocal'*temp8*mphitum! *mphi ' *temp8*mphi; 


EETECEVEEEEEEEEEEETELEE EE EEEETEEETEBEEETEEETELEEEEEEEEETELEEEEEEEEEEEEEEEEEEES 


% % 
% This function creates the product of matrices for integration of FQ % 
% coefficient matrix. % 
% % 


EEGCEEEEEEEFEEEEEEEEEEEBEEETEEEEEBEEEEETEEEEEEELETETELETEEEEEE TETEBETELETETETS 


function y33=astring33 (x,c1,c2,betal, beta2, f1, £2, £3, f£4,mphi,w,... 
wresid, rlocal, templ1, temp2, um, tenp, rlocal, wdot, umdot, gravvec, temp7,... 
tenp8, temp9, temp10, mobhi2, stif) 


mohi (3, 1)=f£1* (c1* (cos (betal*x) +cosh (betal*x) )+sin (betal*x)+... 
Sinh (betal*x) )+£3* (c2* (cos (beta2*x) +cosh (beta2*x))+... 
sin (beta2*x) +sinh (beta2*x) ); 


mohi (3, 2) =£2* (c1* (cos (betal *x) +cosh (betal*x))+sin(betal*x)+... 
sinh (betal*x) ) +£4* (c2* (cos (beta2*x) +cosh (beta2*x))+... y 
sin (beta2*x) +sinh (beta2*x) ); 


rlocal (1) =1.0; 
rlocal (2) =x; 


y33=-rlocal '*temp9* rlocal-rilocal ' *temp9*mphi *um-2* rlocal '*temp10*... 
mpl i *umdot—um! *mphi ' *temp9* rlocal-um' *mphi ' *temp9*nphi*um-2*... 

um’ *mphi ' *temp10*mphi*umdot+rlocal'*temp' *gravvectgravvec' *temp*... 
moh *um; 


EETCTETEEEEEEETETEEEEEEEETEEETEETETETETETETETEEEEEEEEELELE LE TELE EEEEEEEEEEEE SES 
% % 
% This function creates the product of matrices for integration of MN % 
% coefficient matrix. % 
% % 
TELETTEEEETELEEEETEEEEEEEEEELELEEEEESESEEEBELEEEETTELEEEEEEEBETEEEEEEEETEEEEEEES 


function y44=astring44 (x, cl, c2,betal, beta2, f1, £2, £3, £4,mphi,w,... 
wresid, rlocal, templ, temp2, um, temp, rlocal, wdot, umdot, gravvec, temp7,... 
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temp8, temp9, temp10, mphi2, stif) 


mphi (3, 1) =£1* (cl* (cos (betal*x) +cosh (betal*x))4sin(betal*x)+... 
sinh (betal*x) ) +£3* (c2* (cos (beta2*x) tcosh (beta2*x)) +... 


sin (beta2*x) +sinh (beta2*x) ); 


inphi (3, 2) =£2* (cl* (cos (betal*x) tcosh (betal*x))+sin(betal*x)+... 
sinh (betal*x) )+£4* (c2* (cos (beta2*x) +cosh (beta2*x))+... 
sin (beta2*x) +sinh (beta2*x) ); 


v44=mphi !*w'! *w*mphi; 
CEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEELEEEEEEEELEEEEEEELELEBEEE ES 
3 % 
% This function creates the product of matrices for integration of GN % 
% coefficient matrix. % 
% 


ol? 


SEEEEEEEEEELEEELEEEELEEEEEEEEEEEEEEEEEEEEEEEEEEEEEELEEELELELEEEEEEEEEEELE TEE 


ge 


function y55=astring55 (x,cl,c2,betal, beta2, f1, £2, £3, £4,mphi,w,... 
wresid, rlocal, templ1, tenp2, um, tenp, rlocal, wdot, umdot, gravvec, temp7,... 


temp8, temp9, templ0, mphi2, stif) 


mphi (3, 1)=f£1* (cl* (cos (betal*x) +cosh (betal*x) )+sin (betal*x)+... 
sinh (betal*x) ) +£3* (c2* (cos (beta2*x) +cosh (beta2*x))+... 


_ sin (beta2*x) +sinh (beta2*x) ); 


rphi (3, 2) =£2* (cl* (cos (betal*x) +cosh (betal*x) )+sin (betal*x) +... 
sinh (betal*x) )+£4* (c2* (cos (beta2*x) +cosh (beta2*x)) +... y 
sin (beta2*x) +sinh (beta2*x) ); 


yoS=2*mphi' *w' *wdot*mphi; 


% 


% This function creates the product of matrices for integration of KN % 
% coefficient matrix. % 
% 


ESEEEEEEEEEEEEETEEELELEEEEEEEBEEEEEEEE EEE EEE GEEEEEEEEEEEEEEEEEEEEEEEEEELEEE EES 


% 


function y66=astring66 (x,cl,c2,betal, beta2, f1, £2, £3, £4,mphi,w,... 
wresid, rlocal, temp1, temp2, um, tenp, rlocal, wdot, umdot, gravvec, temp7,... 


tarps, temp9, temp10, mphi2, stif) 
moti (3, 1) =£1* (cl* (cos (betal*x) +cosh (betal*x) ) tsin (betal*x) +... 


sinh (betal*x) ) +£3* (c2* (cos (beta2*x) +cosh (beta2*x))+... 
Sin (beta2*x) +sinh (beta2*x) ); 
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mplii (3,2) =£2* (c1* (cos (betal*x) +cosh (betal*x) )+sin (betal*x)+... 
sinh (betal*x) )+£4* (c2* (cos (bet a2*x) +cosh (beta2*x))+... 

sin (beta2*x) +sinh (beta2*x) )e- 

y66-mphi '*w' *wresid*mphi; 


EEEUTEBTEEEELEEEBETEEETEEETE TEBE TE TEE EETEEEEEEEEETEEEEETEEEEEEEEEE EEEEEEEEEEE EES 


$ % 
% This function creates the product of matrices for integration of KN % 
% coefficient matrix. % 
% % 


EEG SSEEEEEEEEEEEEEETEEEEEEEEEEEEEEEEEEEEEEEEEEEEEETEEEEEEEEEEEEEEEEEEEEEEEEEEES 


function y77=astring?7 (x,cl,c2,betal, beta2, f1, £2, £3, £4,mphi,w,... 
wresid, rlocal, temp1, temp2, um, temp, rlocal, wdot, umdot, gravvec, temp?7,... 
tenp8, temp9, temp10, mphi2, stif) 


mohi2 (3, 1)=£1* (betal%2) * (c1* (-cos (betal*x) +cosh (betal*x))... 
~sin (betal*x) +sinh (betal*x) )+£3* (beta2%2) * (c2*... 
(-cos (beta2*x) +cosh (beta2*x) ) -sin (beta2*x) +sinh (beta2*x) ) ; 


moli2 (3,2) =£2* (betal%*2) * (c1l* (-cos (beta1*x) 4cosh (betal*x))... 
-sin (betal*x) +sinh (betal*x) )+£4* (beta2%2) * (c2*... 
(~cos (beta2*x) +cosh (beta2*x) )-sin (beta2*x) +sinh (beta2*x) ); 


y/7J=nphi2'*stif*mphi2; 

ESVETEREBEEEBEBSEEEEEEBRBRBRBBEEELBEREREEBSEEEIFETISEETESSEBEESEEBSEEIGESSBSI 
% J % 
% This function creates the product of matrices for integration of FN % 
% coefficient matrix. % 
% % 


SEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEELEEEEELEEEEEEELELEEELELEELEEESE 


function y88=astring88 (x,cl,c2,betal, beta2, f1,f2, £3, £4,mphi,w,... 
wresid, rlocal, temp1, temp2, um, temp, rlocal, wdot, umdot, gravvec,temp7,... 
tenp8, temp9, temp10, mphi2, stif) 

mohi (3,1) =£1* (cl* (cos (betal*x) +cosh (betal*x) ) +sin (betal*x)+... 

sinh (betal*x) ) +£3* (c2* (cos (beta2*x) +cosh (beta2*x)) +... 

sin (beta2*x) +sinh (beta2*x) ) ; 

mpohi (3,2) =£2* (c1* (cos (betal*x) +cosh (betal*x) ) tsin (betal*x) +... 

sinh (betal*x) ) +£4* (c2* (cos (beta2*x) +cosh (beta2*x))+... 

Sin (beta2*x) +sinh (beta2*x) ) ; 


rlocal(1)=1; 
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rlocal (2) =x; 
5 


y88=mphi'*w' *gravvec-impohi'*w' *wresid*rlocal; 
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